This logbook will answer the Exercises of the module “Tree phenology analysis in R”. It will start with chapter 3, ends with chapter 35 and will only shortly tell if a chapter is left out (like chapter 17).

Chapter 03 Tree dormancy

  1. Put yourself in the place of a breeder who wants to calculate the temperature requirements of a newly released cultivar. Which method will you use to calculate the chilling and forcing periods?

    • Because the cultivar is already released there should be official data about the cultivar. Therefore I would use these data to calculate the temperature requirements. As base data I need the long phenological data sets and temperature records. With these I can calculate finally the date of dormancy overcome.
  2. Which are the advantages (name 2) of the BBCH scale compared with early scales?

    • Growth stages are easily recognizable under field conditions.

    • Growth stages are graded in the order of appearance.

  3. Classify the following phenological stages of sweet cherry according to the BBCH scale.

    The pictures are described from left to right. All of the classifications are based on the given pictures and I assume that the image section shown corresponds to the most represented of the tree.

    • Picture 1: BBCH-stage 03
      • Main stage: 00: sprouting
      • Stage 03: end of the bud swelling (leave buds), parts of the buds are light green, the bud scales are spread apart
    • Picture 2: BBCH-stage 65
      • Main stage 60: Flowering
      • Stage 65: Full flowering, at least 50 % of the flowers are open, the first petals are already falling.
    • Picture 3: BBCH-stage: 86
      • Main stage 80: ripening or maturity
      • Stage 85: Coloring advanced (The color is too light to be ripe but is clearly advanced).

phenological stages of cherry Picture 01: Phenological stages of a cherry tree

Chapter 04 Climate change and impact projection

  1. List the main drivers of climate change at the decade to century scale, and briefly explain the mechanism through which the currently most important driver affects our climate.
    • Main drivers of climate change (influence from decade to century scale)
      • Greenhouse gases
      • Aerosols
      • Clouds
      • Surface albedo
      • Ozone
      • Sun
      • Long-term drivers
    • Currently most important for the climate are greenhouse gases and their increasing amount in the atmosphere
      • The reason are mainly human activities like deforestation, fossil fuel combustion and industrial processes

      • Greenhouse gases trap heat in the atmosphere

      • Can absorb only radiation of certain wavelengths

      • Short-wave radiation from Earth surface is absorbed

  2. Explain briefly what is special about temperature dynamics of recent decades, and why we have good reasons to be concerned.
    • Recently the mean temperatures are globally rising clearly and very drastically

    • Especially in the last 50 to 100 years

    • This an other factors have a great influence in our climate and weather (regional and global)

    • We will and already do experience more extreme weather (like more dry periods, extreme rain, especially mild and extreme winter)

    • Most important is the rapid warming (as conclusion)

  3. What does the abbreviation ‘RCP’ stand for, how are RCPs defined, and what is their role in projecting future climates?
    • RCP -> Representative Concentration Pathways

    • RCPs project future greenhouse gas concentrations (NOT emissions) based on different climate change scenarios until the year 2100

    • The different scenarios are labeled after the possible radioative forcing values in the year 2100

  4. Briefly describe the 4 climate impact projection methods described in the fourth video.
    • Statistical models
      • Analyzes relationships between species occurrences and environmental factors to predict distributions. Often based on historical data.
      • Does not account for ecological interactions or future adaptability.
    • Species distribution modeling
      • Uses species presence/absence data and environmental variables to estimate potential distributions under current and future climates.
      • Assumes species are in equilibrium with their environment, which may hold under climate change.
    • Process-based models.
      • Simulate biological and ecological processes (e.g. metabolism, reproduction) to predict species and ecosystem responses.
      • Requires extensive species-specific data and is computationally demanding.
    • Climate analogue analysis
      • Identifies regions where current conditions resemble projected future climates to predict potential species shifts.

      • Does not consider non-climatic factors (e.g. land use, human impact) that influence species survival.

Chapter 05 Winter chill projections

  1. Sketch out three data access and processing challenges that had to be overcome in order to produce chill projections with state-of-the-art methodology.
    • Climate Data Acquisition
      • Climate datasets come from multiple sources, such as meteorological stations, reanalysis datasets and climate models.
      • Differences in spatial and temporal resolution require preprocessing to ensure consistency.
      • Missing data and measurement errors need handling to avoid distortions in projections.
    • Model Selection and Calibration
      • Multiple winter chill models exist (e.g. Chill Hours, Utah Model, Dynamic Model) each responding differently to temperature variations.
      • Choosing the most appropriate model for a specific region is crucial, as different tree species have varying chill requirements.
      • Calibration is needed to align models with observed historical data and improve predictive accuracy.
    • Uncertainty and Bias Correction
      • Climate models often introduce biases that must be corrected before using their projections.

      • Different greenhouse gas emission scenarios (e.g. RCP2.6, RCP4.5, RCP8.5) result in a range of possible future conditions, adding uncertainty.

      • Statistical methods, such as ensemble modeling and probabilistic approaches, help quantify and communicate these uncertainties.

  2. Outline, in your understanding, the basic steps that are necessary to make such projections.
    1. Data collection and Preprocessing

      • Gather historical temperature records and future climate projections from global and regional sources.
      • Standardize the data sets by converting them into compatible formats and resolutions.
      • Fill in missing data and remove inconsistencies that could affect model accuracy.
    2. Chill Model Application

      • Apply a chosen winter chill model to calculate chill accumulation over time.
      • Compare results from different models to understand their variability and sensitivity.
      • Validate model outputs against historical observations to ensure reliability.
    3. Climate Scenario Analysis

      • Run chill models using projected temperature data under different climate change scenarios.
      • Analyse how winter chill accumulation might shift in various regions over time.
      • Compare results across different periods (e.g. near-future vs. late-century) to assess long-term trends.
    4. Validation and Uncertainty Assessment

      • Compare modeled projections with historical observations to evaluate accuracy.

      • Apply statistical techniques to quantify uncertainties and assess model reliability.

      • Use ensemble approaches to minimize errors and account for variability between models.

    5. Visualization and Interpretation

      • Present results using graphs, maps, and reports for easy interpretation.

      • Provide insights for stakeholders, such as farmers and policymakers, to support decision-making.

      • Communicate potential risks and adaptation strategies based on projected changes.

Chapter 06 Manual chill analysis

  1. Write a basic function that calculates warm hours (>25° C).

    • Basic structure:

      • # read in Dataset (e.g. Winters_hours_gaps)

      • dataset[ , new_column_name] <- dataset$temperature_column > 25

  2. Apply this function to the Winters_hours_gaps dataset.

Winters_hours_gaps <- Winters_hours_gaps
hourtemps <- Winters_hours_gaps[ , c("Year", "Month", "Day", "Hour", "Temp")]

hourtemps[ , "warm_hour"] <- hourtemps$Temp > 25

hourtemps[2503:2507, ]
##      Year Month Day Hour   Temp warm_hour
## 2503 2008     6  15   16 27.604      TRUE
## 2504 2008     6  15   17 26.989      TRUE
## 2505 2008     6  15   18 24.436     FALSE
## 2506 2008     6  15   19 23.088     FALSE
## 2507 2008     6  15   20 19.555     FALSE
  1. Extend this function, so that it can take start and end dates as inputs and sums up warm hours between these dates.
Start_Date <- which(hourtemps$Year == 2008 &
                      hourtemps$Month == 6 &
                      hourtemps$Day == 1 &
                      hourtemps$Hour == 12)

End_date <- which(hourtemps$Year == 2008 &
                    hourtemps$Month == 6 &
                    hourtemps$Day == 30 &
                    hourtemps$Hour ==12)
sum(hourtemps$warm_hour[Start_Date:End_date])
## [1] 250
# The result is 250 and will be the control for the following function with the same start and end

sum_warmH <- function(WHourly, 
                        startYear,
                        startMonth,
                        startDay,                   
                        startHour,
                        endYear,
                        endMonth,
                        endDay,                   
                        endHour)
{WHourly[,"warm_hour2"] <- WHourly$Temp > 25
  Start_row <- which(hourtemps$Year == startYear &
                     hourtemps$Month == startMonth &
                     hourtemps$Day == startDay &
                     hourtemps$Hour == startHour)

  End_row <- which(hourtemps$Year == endYear &
                   hourtemps$Month == endMonth &
                   hourtemps$Day == endDay &
                   hourtemps$Hour == endHour)

WHs <- sum(WHourly$warm_hour[Start_row:End_row])
    
return(WHs)}

sum_warmH(hourtemps,
           startYear = 2008,
           startMonth = 6,
           startDay = 1,
           startHour = 12,
           endYear = 2008,
           endMonth = 6,
           endDay = 30,
           endHour = 12)
## [1] 250

Chapter 07 Chill models

  1. Run the “chilling()” function on the Winters_hours_gap dataset.
hourtemps <- Winters_hours_gaps[ , c("Year", "Month", "Day", "Hour", "Temp")]

chilling_l7 <- chilling(make_JDay(hourtemps),
                   Start_JDay = 90,
                   End_JDay = 100)
chilling_l7
##      Season End_year Season_days Data_days Perc_complete Chilling_Hours
## 1 2007/2008     2008          11        11           100             40
##   Utah_Model Chill_portions     GDH
## 1       15.5       2.009147 2406.52
  1. Create your own temperature-weighting chill model using the “step_model()” function.
df_l7 = data.frame(
  lower = c(-1000, 1, 2, 3, 4, 5, 6),
  upper = c(    1, 2, 3, 4, 5, 6, 1000),
  weight = c(   0, 1, 2, 3, 2, 1, 0))
  1. Run this model on the Winters_hours_gap dataset using the “tempResponse()” function.
step_model_l7 <- step_model(HourTemp=hourtemps$Temp, df=df_l7)

TResp_l7 <- tempResponse(make_JDay(hourtemps),
                         Start_JDay = 90,
                         End_JDay = 100,
                         models = list(Chill_Portions = Dynamic_Model,
                                       GDH = GDH))
TResp_l7
##      Season End_year Season_days Data_days Perc_complete Chill_Portions     GDH
## 1 2007/2008     2008          11        11           100       2.009147 2406.52

Chapter 08 Making hourly temperatures

  1. Choose a location of interest, find out its latitude and produce plots of daily sunrise, sunset and day length.
    • Location of interest: Dairy farm with own fodder cultivation (belongs to my uncle)
    • Bislich, 46487 Wesel
    • Latitude: 51°40’58.4”N Longitude: 6°29’46.4”E
Days_l8 <- daylength(latitude = 51.4, JDay = 1:365)
Days_df_l8 <- data.frame(JDay = 1:365,
                         Sunrise = Days_l8$Sunrise,
                         Sunset = Days_l8$Sunset,
                         Daylength = Days_l8$Daylength)

Days_df_l8 <- pivot_longer(Days_df_l8, cols=c(Sunrise:Daylength))

ggplot(Days_df_l8, aes(JDay, value)) +
        geom_line(lwd = 1.5) +
        facet_grid(cols = vars(name)) +
        ylab("Time of Day / Daylength [hours]") +
        xlab("Day of the year [JDay]") +
        theme_bw(base_size = 20)

  1. Produce an hourly dataset, based on idealized daily curves, for the KA_weather dataset (included in ChillR).
hourtemp_KA <- stack_hourly_temps(KA_weather, latitude=50.4)
  1. Produce empirical temperature curve parameters for the Winters_hours_gaps dataset, and use them to predict hourly values from daily temperatures (this is very similar to the example above, but please make sure you understand what’s going on).
empi_curve_l8 <- Empirical_daily_temperature_curve(Winters_hours_gaps)

ggplot(data = empi_curve_l8[1:96, ], aes(Hour, Prediction_coefficient)) +
  geom_line(lwd = 1.3, col = "red") +
  facet_grid(rows = vars(Month)) +
  xlab("Hour of the day") +
  ylab("Prediction coefficient") +
  theme_bw(base_size = 20)

coeffs_l8 <- Empirical_daily_temperature_curve(Winters_hours_gaps)
Winters_daily_l8 <- make_all_day_table(Winters_hours_gaps, input_timestep = "hour")
Winter_hours_l8 <- Empirical_hourly_temperatures(Winters_daily_l8, coeffs_l8)


# simplify the data for easier use
Winter_hours_l8 <- Winter_hours_l8[ , c("Year", "Month", "Day", "Hour", "Temp")]
colnames(Winter_hours_l8)[ncol(Winter_hours_l8)] <- "Temp_empirical"
Winters_ideal_l8 <- stack_hourly_temps(Winters_daily_l8, latitude = 38.5)$hourtemps
Winters_ideal_l8 <- Winters_ideal_l8[ , c("Year", "Month", "Day", "Hour", "Temp")]
colnames(Winters_ideal_l8)[ncol(Winters_ideal_l8)] <- "Temp_ideal"


# create the "triangular" data set to compare in the end
Winters_triangle_l8 <- Winters_daily_l8
Winters_triangle_l8[ , "Hour"] <- 0
Winters_triangle_l8$Hour[nrow(Winters_triangle_l8)] <- 23
Winters_triangle_l8[ , "Temp"] <- 0
Winters_triangle_l8 <- make_all_day_table(Winters_triangle_l8, timestep = "hour")
colnames(Winters_triangle_l8)[ncol(Winters_triangle_l8)] <- "Temp_triangular"


# the following loop will fill in the daily Tmin and Tmax values for every hour
for (i in 2:nrow(Winters_triangle_l8))
{if (is.na(Winters_triangle_l8$Tmin[i]))
  Winters_triangle_l8$Tmin[i] <- Winters_triangle_l8$Tmin[i - 1]
 if (is.na(Winters_triangle_l8$Tmax[i]))
   Winters_triangle_l8$Tmax[i] <- Winters_triangle_l8$Tmax[i - 1]}

Winters_triangle_l8$Temp_triangular <- NA


# Tmin is fixated to the 6th hour of each day
Winters_triangle_l8$Temp_triangular[which(Winters_triangle_l8$Hour == 6)] <-
  Winters_triangle_l8$Tmin[which(Winters_triangle_l8$Hour == 6)]


# Tmax is fixated to the 18th hour of each day
Winters_triangle_l8$Temp_triangular[which(Winters_triangle_l8$Hour == 18)] <-
  Winters_triangle_l8$Tmax[which(Winters_triangle_l8$Hour == 18)]


# the following loop will fill in the gaps between the min and max data in a linear way
Winters_triangle_l8$Temp_triangular <- interpolate_gaps(Winters_triangle_l8$Temp_triangular)$interp
Winters_triangle_l8 <- Winters_triangle_l8[ , c("Year", "Month", "Day", "Hour", "Temp_triangular")]


# merge all created data frames for easy display
Winters_temps_l8 <- merge(Winters_hours_gaps, Winter_hours_l8,
                          by = c("Year", "Month", "Day", "Hour"))
Winters_temps_l8 <- merge(Winters_temps_l8, Winters_triangle_l8,
                          by = c("Year", "Month", "Day", "Hour"))
Winters_temps_l8 <- merge(Winters_temps_l8, Winters_ideal_l8,
                          by = c("Year", "Month", "Day", "Hour"))


# convert the date into R's date format and reorganize the data frame
Winters_temps_l8[, "DATE"] <-
  ISOdate(Winters_temps_l8$Year,
          Winters_temps_l8$Month,
          Winters_temps_l8$Day,
          Winters_temps_l8$Hour)


Winters_temps_to_plot_l8 <-
  Winters_temps_l8[, c("DATE",
                    "Temp",
                    "Temp_empirical",
                    "Temp_triangular",
                    "Temp_ideal")]
Winters_temps_to_plot_l8 <- Winters_temps_to_plot_l8[100:200, ]
Winters_temps_to_plot_l8 <- pivot_longer(Winters_temps_to_plot_l8, cols=Temp:Temp_ideal)
colnames(Winters_temps_to_plot_l8) <- c("DATE", "Method", "Temperature")


ggplot(data = Winters_temps_to_plot_l8, aes(DATE, Temperature, colour = Method)) +
  geom_line(lwd = 1.3) + ylab("Temperature (°C)") + xlab("Date")

Chapter 09 Some useful tools in R

  1. Based on the Winters_hours_gaps data set, use magrittr pipes and functions of the tidyverse to accomplish the following:
    1. Convert the data set into a tibble
    2. Select only the top 10 rows of the data set
    3. Convert the tibble to a long format, with separate rows for Temp_gaps and Temp
    4. Use ggplot2 to plot Temp_gaps and temp as facets (point or line plot)
    5. Convert the data set back to the wide format
    6. Select only the following columns: Year, Month, Day and Temp
    7. Sort the data set by the Temp column, in descending order
# Convert the data set into a tibble
hourtemps <- Winters_hours_gaps
hourtemps %>% as_tibble() %>% summary() 
##       Year          Month             Day             Hour     
##  Min.   :2008   Min.   : 3.000   Min.   : 1.00   Min.   : 0.0  
##  1st Qu.:2008   1st Qu.: 5.000   1st Qu.: 8.00   1st Qu.: 6.0  
##  Median :2008   Median : 7.000   Median :15.00   Median :11.0  
##  Mean   :2008   Mean   : 6.722   Mean   :15.53   Mean   :11.5  
##  3rd Qu.:2008   3rd Qu.: 9.000   3rd Qu.:23.00   3rd Qu.:17.0  
##  Max.   :2008   Max.   :11.000   Max.   :31.00   Max.   :23.0  
##                                                                
##    Temp_gaps           Temp       
##  Min.   : 1.967   Min.   : 1.967  
##  1st Qu.:13.016   1st Qu.:13.257  
##  Median :17.915   Median :18.010  
##  Mean   :18.419   Mean   :18.644  
##  3rd Qu.:23.569   3rd Qu.:23.857  
##  Max.   :38.170   Max.   :38.365  
##  NA's   :3374
# Select only the top 10 rows of the data set
hourtemps_l9 <- as_tibble(hourtemps[1:10, ])

# Convert the tibble to a long format, with separate rows for Temp_gaps and Temp
hourtemps_long_l9 <- hourtemps_l9 %>% pivot_longer(cols = Temp:Temp_gaps)


# Use ggplot2 to plot Temp_gaps and temp as facets (point or line plot)
ggplot(hourtemps_l9, aes(Hour, Temp_gaps)) + geom_point()

# Convert the data set back to the wide format
hourtemps_wide_l9 <- hourtemps_long_l9 %>% pivot_wider(names_from = name)

# Select only the following columns: Year, Month, Day and Temp
hourtemps_wide_l9 %>% select(c(Month, Day, Temp))
## # A tibble: 10 × 3
##    Month   Day  Temp
##    <int> <int> <dbl>
##  1     3     3  15.1
##  2     3     3  17.2
##  3     3     3  18.7
##  4     3     3  18.7
##  5     3     3  18.8
##  6     3     3  19.5
##  7     3     3  19.3
##  8     3     3  17.7
##  9     3     3  15.4
## 10     3     3  12.7
# Sort the data set by the Temp column, in descending order
hourtemps_wide_l9 %>% arrange(desc(Temp))
## # A tibble: 10 × 6
##     Year Month   Day  Hour  Temp Temp_gaps
##    <int> <int> <int> <int> <dbl>     <dbl>
##  1  2008     3     3    15  19.5      19.5
##  2  2008     3     3    16  19.3      19.3
##  3  2008     3     3    14  18.8      18.8
##  4  2008     3     3    12  18.7      18.7
##  5  2008     3     3    13  18.7      18.7
##  6  2008     3     3    17  17.7      17.7
##  7  2008     3     3    11  17.2      17.2
##  8  2008     3     3    18  15.4      15.4
##  9  2008     3     3    10  15.1      15.1
## 10  2008     3     3    19  12.7      12.7
  1. For the Winter_hours_gaps data set, write a for loop to convert all temperatures (Temp column) to degrees Fahrenheit.
HT_l9.2 <- Winters_hours_gaps
HT_l9.2$Temp_F <- NA
for (i in 1:nrow(HT_l9.2))
  {HT_l9.2$Temp_F[i] <- (HT_l9.2$Temp[i] * 9/5) + 32}
head(HT_l9.2)
##   Year Month Day Hour Temp_gaps   Temp  Temp_F
## 1 2008     3   3   10    15.127 15.127 59.2286
## 2 2008     3   3   11    17.153 17.153 62.8754
## 3 2008     3   3   12    18.699 18.699 65.6582
## 4 2008     3   3   13    18.699 18.699 65.6582
## 5 2008     3   3   14    18.842 18.842 65.9156
## 6 2008     3   3   15    19.508 19.508 67.1144
  1. Execute the same operation with a function from the apply family.
HT_l9.3 <- hourtemps
func_Far_l9.3 <- function(x)
  {(x * 9/5) + 32}
head(HT_l9.3)
##   Year Month Day Hour Temp_gaps   Temp
## 1 2008     3   3   10    15.127 15.127
## 2 2008     3   3   11    17.153 17.153
## 3 2008     3   3   12    18.699 18.699
## 4 2008     3   3   13    18.699 18.699
## 5 2008     3   3   14    18.842 18.842
## 6 2008     3   3   15    19.508 19.508
for (i in 1:nrow(hourtemps))
  {HT_l9.3$Temp_F[i] <- sapply(HT_l9.3$Temp[i], func_Far_l9.3)}
head(HT_l9.3)
##   Year Month Day Hour Temp_gaps   Temp  Temp_F
## 1 2008     3   3   10    15.127 15.127 59.2286
## 2 2008     3   3   11    17.153 17.153 62.8754
## 3 2008     3   3   12    18.699 18.699 65.6582
## 4 2008     3   3   13    18.699 18.699 65.6582
## 5 2008     3   3   14    18.842 18.842 65.9156
## 6 2008     3   3   15    19.508 19.508 67.1144
  1. Now use the tidyverse function mutate to achieve the same outcome
HT_l9.2 <- hourtemps %>% mutate(Temp_F = (Temp * 9/5) + 32)
head(HT_l9.2)
##   Year Month Day Hour Temp_gaps   Temp  Temp_F
## 1 2008     3   3   10    15.127 15.127 59.2286
## 2 2008     3   3   11    17.153 17.153 62.8754
## 3 2008     3   3   12    18.699 18.699 65.6582
## 4 2008     3   3   13    18.699 18.699 65.6582
## 5 2008     3   3   14    18.842 18.842 65.9156
## 6 2008     3   3   15    19.508 19.508 67.1144

Chapter 10 Getting temperature data

  1. Choose a location of interest and find the 25 closest weather stations using the handle_gsod function.
  • location of interest: Dairy farm with own fodder cultivation (belongs to my uncle) Bislich, 46487 Wesel
  • coordinates: latitude: 51°40’58.4”N longitude: 6°29’46.4”E
station_list <- handle_gsod(action = "list_stations",
                            location = c(6.29, 51.40),
                            time_interval = c(1990, 2020))
station_list <- slice(station_list, 1:25)

head(station_list)
## # A tibble: 6 × 10
##   chillR_code STATION.NAME        CTRY    Lat  Long    BEGIN      END Distance
##   <chr>       <chr>               <chr> <dbl> <dbl>    <int>    <int>    <dbl>
## 1 06391099999 ARCEN AWS           "NL"   51.5  6.2  19940901 20241127     12.8
## 2 10403099999 MOENCHENGLADBACH    "GM"   51.2  6.5  19381001 19421031     23.6
## 3 10437499999 MONCHENGLADBACH     "GM"   51.2  6.50 19960715 20241127     24.1
## 4 10405099999 LAARBRUCH (RAF)     "GM"   51.6  6.15 19730102 19971231     24.3
## 5 10000199999 NIEDERRHEIN         "GM"   51.6  6.14 20050101 20241127     24.7
## 6 99860099999 ARRC MOENCHENGLADBA ""     51.2  6.21 20011015 20020306     24.7
## # ℹ 2 more variables: Overlap_years <dbl>, Perc_interval_covered <dbl>
  1. Download weather data for the most promissing station on the list.
weather_WES <- handle_gsod(action = "download_weather",
                          location = station_list$chillR_code[1],
                          time_interval = c(1990, 2020))
  1. Convert the weather data into chillR format and save the data on your PC.
cleaned_weather_WES <- handle_gsod(weather_WES)
cleaned_weather_WES[[1]][1000:1010,]
write.csv(station_list, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/station_list_Wesel.csv", row.names = FALSE)
write.csv(weather_WES[[1]], "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_raw_weather.csv", row.names = FALSE)
write.csv(cleaned_weather_WES[[1]], "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_chillR_weather.csv", row.names = FALSE)

Chapter 11 Filling gaps in temperature records

Before the answers of the exercises for this chapter start, there will be dates added to the dataset because the first 4 years are missing for a reason. So the get added here manually to answer all following exercises correctly.

cleaned_weather_WES <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_chillR_weather.csv")

start_date <- as.Date("1990-01-01")
end_date <- as.Date("1993-12-31")
date_sequence <- seq.Date(start_date, end_date, by = "day")

new_rows <- data.frame(
  Date = as.POSIXct(paste(date_sequence, "12:00:00")),
  Year = as.numeric(format(date_sequence, "%Y")),
  Month = as.numeric(format(date_sequence, "%m")),
  Day = as.numeric(format(date_sequence, "%d")),
  Tmin = NA,
  Tmax = NA,
  Tmean = NA,
  Prec = NA
)
cleaned_weather_WES <- rbind(new_rows, cleaned_weather_WES)
  1. Use chillR functions to find out how many gaps you have in this dataset (even if you have none, please still follow all further steps).
Wesel_QC <- fix_weather(cleaned_weather_WES)$QC

table(is.na(cleaned_weather_WES$Tmin))
## 
## FALSE  TRUE 
##  9022  2301
Wesel_QC
##       Season End_year Season_days Data_days Missing_Tmin Missing_Tmax
## 1  1989/1990     1990         365       365          365          365
## 2  1990/1991     1991         365       365          365          365
## 3  1991/1992     1992         366       366          366          366
## 4  1992/1993     1993         365       365          365          365
## 5  1993/1994     1994         365       365          270          270
## 6  1994/1995     1995         365       365          161          161
## 7  1995/1996     1996         366       366          128          128
## 8  1996/1997     1997         365       365          152          152
## 9  1997/1998     1998         365       365           46           46
## 10 1998/1999     1999         365       365            6            6
## 11 1999/2000     2000         366       366            0            0
## 12 2000/2001     2001         365       365           31           31
## 13 2001/2002     2002         365       365            4            4
## 14 2002/2003     2003         365       365            2            2
## 15 2003/2004     2004         366       366            2            2
## 16 2004/2005     2005         365       365           11           11
## 17 2005/2006     2006         365       365            2            2
## 18 2006/2007     2007         365       365            4            4
## 19 2007/2008     2008         366       366            5            5
## 20 2008/2009     2009         365       365            0            0
## 21 2009/2010     2010         365       365            0            0
## 22 2010/2011     2011         365       365            0            0
## 23 2011/2012     2012         366       366            4            4
## 24 2012/2013     2013         365       365            5            5
## 25 2013/2014     2014         365       365            2            2
## 26 2014/2015     2015         365       365            2            2
## 27 2015/2016     2016         366       366            0            0
## 28 2016/2017     2017         365       365            0            0
## 29 2017/2018     2018         365       365            2            2
## 30 2018/2019     2019         365       365            0            0
## 31 2019/2020     2020         366       366            1            1
##    Incomplete_days Perc_complete
## 1              365           0.0
## 2              365           0.0
## 3              366           0.0
## 4              365           0.0
## 5              270          26.0
## 6              161          55.9
## 7              128          65.0
## 8              152          58.4
## 9               46          87.4
## 10               6          98.4
## 11               0         100.0
## 12              31          91.5
## 13               4          98.9
## 14               2          99.5
## 15               2          99.5
## 16              11          97.0
## 17               2          99.5
## 18               4          98.9
## 19               5          98.6
## 20               0         100.0
## 21               0         100.0
## 22               0         100.0
## 23               4          98.9
## 24               5          98.6
## 25               2          99.5
## 26               2          99.5
## 27               0         100.0
## 28               0         100.0
## 29               2          99.5
## 30               0         100.0
## 31               1          99.7
  1. Create a list of the 25 closest weather stations using the handle_gsod function.
station_list <- handle_gsod(action = "list_stations",
                            location = c(6.29, 51.40),
                            time_interval = c(1990, 2020))
station_list <- slice(station_list, 1:25)
station_list
## # A tibble: 25 × 10
##    chillR_code STATION.NAME        CTRY    Lat  Long    BEGIN      END Distance
##    <chr>       <chr>               <chr> <dbl> <dbl>    <int>    <int>    <dbl>
##  1 06391099999 ARCEN AWS           "NL"   51.5  6.2  19940901 20241127     12.8
##  2 10403099999 MOENCHENGLADBACH    "GM"   51.2  6.5  19381001 19421031     23.6
##  3 10437499999 MONCHENGLADBACH     "GM"   51.2  6.50 19960715 20241127     24.1
##  4 10405099999 LAARBRUCH (RAF)     "GM"   51.6  6.15 19730102 19971231     24.3
##  5 10000199999 NIEDERRHEIN         "GM"   51.6  6.14 20050101 20241127     24.7
##  6 99860099999 ARRC MOENCHENGLADBA ""     51.2  6.21 20011015 20020306     24.7
##  7 10401099999 BRUGGEN (RAF)       "GM"   51.2  6.13 19730102 19971231     24.8
##  8 06385099999 DE PEEL(NAFB)       "NL"   51.6  5.93 19730402 19940611     29.9
##  9 10402099999 WILDENRATH(GAFB)    "GM"   51.1  6.22 19730101 20030508     31.9
## 10 10404399999 KALKAR (MIL COMM)   "GM"   51.7  6.17 19820927 19870304     32.6
## # ℹ 15 more rows
## # ℹ 2 more variables: Overlap_years <dbl>, Perc_interval_covered <dbl>
  1. Identify suitable weather stations for patching gaps.

The stations 3, 5, 11 are suitable.

  1. Download weather data for promising stations, convert them to ChillR format and compile them in a list
patch_weather <- handle_gsod(action = "download_weather",
                             location = as.character(station_list$chillR_code[c(3, 5, 11)]),
                             time_interval = c(1990,2020)) %>% handle_gsod()
## Loading data for 31 years from station 'NIEDERRHEIN'
## ================================================================================
## 
## Loading data for 31 years from station 'DUSSELDORF'
## ================================================================================
## 
## Loading data for 31 years from station 'MONCHENGLADBACH'
## ================================================================================
  1. Use the patch_daily_temperatures function to fill gaps
patched <- patch_daily_temperatures(weather = cleaned_weather_WES,
                                    patch_weather = patch_weather,
                                    max_mean_bias = 1,
                                    max_stdev_bias = 2)
  1. Investigate the results - have all gaps been filled?
patched$statistics[[1]]
##      mean_bias stdev_bias filled gaps_remain
## Tmin    -0.703      1.151     36        2265
## Tmax     0.251      0.731     36        2265
patched$statistics[[2]]
##      mean_bias stdev_bias filled gaps_remain
## Tmin    -0.541      1.293   1098        1167
## Tmax     0.078      1.036   1098        1167
patched$statistics[[3]]
##      mean_bias stdev_bias filled gaps_remain
## Tmin    -1.206      1.705      0        1167
## Tmax     0.220      1.015    256         911
  1. If necessary, repeat until you have a dataset you can work with in further analyses.
station_list
## # A tibble: 25 × 10
##    chillR_code STATION.NAME        CTRY    Lat  Long    BEGIN      END Distance
##    <chr>       <chr>               <chr> <dbl> <dbl>    <int>    <int>    <dbl>
##  1 06391099999 ARCEN AWS           "NL"   51.5  6.2  19940901 20241127     12.8
##  2 10403099999 MOENCHENGLADBACH    "GM"   51.2  6.5  19381001 19421031     23.6
##  3 10437499999 MONCHENGLADBACH     "GM"   51.2  6.50 19960715 20241127     24.1
##  4 10405099999 LAARBRUCH (RAF)     "GM"   51.6  6.15 19730102 19971231     24.3
##  5 10000199999 NIEDERRHEIN         "GM"   51.6  6.14 20050101 20241127     24.7
##  6 99860099999 ARRC MOENCHENGLADBA ""     51.2  6.21 20011015 20020306     24.7
##  7 10401099999 BRUGGEN (RAF)       "GM"   51.2  6.13 19730102 19971231     24.8
##  8 06385099999 DE PEEL(NAFB)       "NL"   51.6  5.93 19730402 19940611     29.9
##  9 10402099999 WILDENRATH(GAFB)    "GM"   51.1  6.22 19730101 20030508     31.9
## 10 10404399999 KALKAR (MIL COMM)   "GM"   51.7  6.17 19820927 19870304     32.6
## # ℹ 15 more rows
## # ℹ 2 more variables: Overlap_years <dbl>, Perc_interval_covered <dbl>
patch_weather2 <- handle_gsod(action = "download_weather",
                             location = as.character(station_list$chillR_code[c(3, 4, 5, 7, 9, 11, 12, 13, 15, 16)]),
                             time_interval = c(1990,2020)) %>%
  handle_gsod()
## Loading data for 31 years from station 'VOLKEL'
## ================================================================================
## 
## Loading data for 31 years from station 'ELL AWS'
## ================================================================================
## 
## Loading data for 31 years from station 'NIEDERRHEIN'
## ================================================================================
## 
## Loading data for 31 years from station 'DUSSELDORF'
## ================================================================================
## 
## Loading data for 31 years from station 'BRUGGEN (RAF)'
## ================================================================================
## 
## Loading data for 31 years from station 'WILDENRATH(GAFB)'
## ================================================================================
## 
## Loading data for 31 years from station 'KALKAR (MIL COMM)'
## ================================================================================
## 
## Loading data for 31 years from station 'LAARBRUCH (RAF)'
## ================================================================================
## 
## Loading data for 31 years from station 'ESSEN/MULHEIM'
## ================================================================================
## 
## Loading data for 31 years from station 'MONCHENGLADBACH'
## ================================================================================
patched_2 <- patch_daily_temperatures(weather = cleaned_weather_WES,
                                    patch_weather = patch_weather2,
                                    max_mean_bias = 1,
                                    max_stdev_bias = 2)

patched_2$statistics[[1]]
##      mean_bias stdev_bias filled gaps_remain
## Tmin     0.665      1.179   2286          15
## Tmax     0.165      0.873   2286          15
patched_2$statistics[[6]]
##      mean_bias stdev_bias filled gaps_remain
## Tmin     0.459      1.428      0           4
## Tmax    -0.109      1.217      0           4
patched_2$statistics[[7]]
##      mean_bias stdev_bias filled gaps_remain
## Tmin    -0.302      1.258      3           1
## Tmax    -0.018      1.211      3           1
post_patch_stats_2 <- fix_weather(patched_2)$QC
post_patch_stats_2
##       Season End_year Season_days Data_days Missing_Tmin Missing_Tmax
## 1  1989/1990     1990         365       365            0            0
## 2  1990/1991     1991         365       365            0            0
## 3  1991/1992     1992         366       366            0            0
## 4  1992/1993     1993         365       365            0            0
## 5  1993/1994     1994         365       365            0            0
## 6  1994/1995     1995         365       365            0            0
## 7  1995/1996     1996         366       366            0            0
## 8  1996/1997     1997         365       365            0            0
## 9  1997/1998     1998         365       365            0            0
## 10 1998/1999     1999         365       365            1            1
## 11 1999/2000     2000         366       366            0            0
## 12 2000/2001     2001         365       365            0            0
## 13 2001/2002     2002         365       365            0            0
## 14 2002/2003     2003         365       365            0            0
## 15 2003/2004     2004         366       366            0            0
## 16 2004/2005     2005         365       365            0            0
## 17 2005/2006     2006         365       365            0            0
## 18 2006/2007     2007         365       365            0            0
## 19 2007/2008     2008         366       366            0            0
## 20 2008/2009     2009         365       365            0            0
## 21 2009/2010     2010         365       365            0            0
## 22 2010/2011     2011         365       365            0            0
## 23 2011/2012     2012         366       366            0            0
## 24 2012/2013     2013         365       365            0            0
## 25 2013/2014     2014         365       365            0            0
## 26 2014/2015     2015         365       365            0            0
## 27 2015/2016     2016         366       366            0            0
## 28 2016/2017     2017         365       365            0            0
## 29 2017/2018     2018         365       365            0            0
## 30 2018/2019     2019         365       365            0            0
## 31 2019/2020     2020         366       366            0            0
##    Incomplete_days Perc_complete
## 1                0         100.0
## 2                0         100.0
## 3                0         100.0
## 4                0         100.0
## 5                0         100.0
## 6                0         100.0
## 7                0         100.0
## 8                0         100.0
## 9                0         100.0
## 10               1          99.7
## 11               0         100.0
## 12               0         100.0
## 13               0         100.0
## 14               0         100.0
## 15               0         100.0
## 16               0         100.0
## 17               0         100.0
## 18               0         100.0
## 19               0         100.0
## 20               0         100.0
## 21               0         100.0
## 22               0         100.0
## 23               0         100.0
## 24               0         100.0
## 25               0         100.0
## 26               0         100.0
## 27               0         100.0
## 28               0         100.0
## 29               0         100.0
## 30               0         100.0
## 31               0         100.0
patched_complete <- patched_2$weather %>% make_all_day_table()

Tmin_int <- interpolate_gaps(patched_complete[,"Tmin"])

patched_complete <- patched_complete %>% mutate(Tmin = Tmin_int$interp,
                              Tmin_interpolated = Tmin_int$missing)

Tmax_int <- interpolate_gaps(patched_complete[,"Tmax"])

patched_complete <- patched_complete %>% mutate(Tmax = Tmax_int$interp,
                              Tmax_interpolated = Tmax_int$missing)

fix_weather(patched_complete)$QC
##       Season End_year Season_days Data_days Missing_Tmin Missing_Tmax
## 1  1989/1990     1990         365       365            0            0
## 2  1990/1991     1991         365       365            0            0
## 3  1991/1992     1992         366       366            0            0
## 4  1992/1993     1993         365       365            0            0
## 5  1993/1994     1994         365       365            0            0
## 6  1994/1995     1995         365       365            0            0
## 7  1995/1996     1996         366       366            0            0
## 8  1996/1997     1997         365       365            0            0
## 9  1997/1998     1998         365       365            0            0
## 10 1998/1999     1999         365       365            0            0
## 11 1999/2000     2000         366       366            0            0
## 12 2000/2001     2001         365       365            0            0
## 13 2001/2002     2002         365       365            0            0
## 14 2002/2003     2003         365       365            0            0
## 15 2003/2004     2004         366       366            0            0
## 16 2004/2005     2005         365       365            0            0
## 17 2005/2006     2006         365       365            0            0
## 18 2006/2007     2007         365       365            0            0
## 19 2007/2008     2008         366       366            0            0
## 20 2008/2009     2009         365       365            0            0
## 21 2009/2010     2010         365       365            0            0
## 22 2010/2011     2011         365       365            0            0
## 23 2011/2012     2012         366       366            0            0
## 24 2012/2013     2013         365       365            0            0
## 25 2013/2014     2014         365       365            0            0
## 26 2014/2015     2015         365       365            0            0
## 27 2015/2016     2016         366       366            0            0
## 28 2016/2017     2017         365       365            0            0
## 29 2017/2018     2018         365       365            0            0
## 30 2018/2019     2019         365       365            0            0
## 31 2019/2020     2020         366       366            0            0
##    Incomplete_days Perc_complete
## 1                0           100
## 2                0           100
## 3                0           100
## 4                0           100
## 5                0           100
## 6                0           100
## 7                0           100
## 8                0           100
## 9                0           100
## 10               0           100
## 11               0           100
## 12               0           100
## 13               0           100
## 14               0           100
## 15               0           100
## 16               0           100
## 17               0           100
## 18               0           100
## 19               0           100
## 20               0           100
## 21               0           100
## 22               0           100
## 23               0           100
## 24               0           100
## 25               0           100
## 26               0           100
## 27               0           100
## 28               0           100
## 29               0           100
## 30               0           100
## 31               0           100
fix_QC <- fix_weather(patched_complete)$QC
fix_QC
##       Season End_year Season_days Data_days Missing_Tmin Missing_Tmax
## 1  1989/1990     1990         365       365            0            0
## 2  1990/1991     1991         365       365            0            0
## 3  1991/1992     1992         366       366            0            0
## 4  1992/1993     1993         365       365            0            0
## 5  1993/1994     1994         365       365            0            0
## 6  1994/1995     1995         365       365            0            0
## 7  1995/1996     1996         366       366            0            0
## 8  1996/1997     1997         365       365            0            0
## 9  1997/1998     1998         365       365            0            0
## 10 1998/1999     1999         365       365            0            0
## 11 1999/2000     2000         366       366            0            0
## 12 2000/2001     2001         365       365            0            0
## 13 2001/2002     2002         365       365            0            0
## 14 2002/2003     2003         365       365            0            0
## 15 2003/2004     2004         366       366            0            0
## 16 2004/2005     2005         365       365            0            0
## 17 2005/2006     2006         365       365            0            0
## 18 2006/2007     2007         365       365            0            0
## 19 2007/2008     2008         366       366            0            0
## 20 2008/2009     2009         365       365            0            0
## 21 2009/2010     2010         365       365            0            0
## 22 2010/2011     2011         365       365            0            0
## 23 2011/2012     2012         366       366            0            0
## 24 2012/2013     2013         365       365            0            0
## 25 2013/2014     2014         365       365            0            0
## 26 2014/2015     2015         365       365            0            0
## 27 2015/2016     2016         366       366            0            0
## 28 2016/2017     2017         365       365            0            0
## 29 2017/2018     2018         365       365            0            0
## 30 2018/2019     2019         365       365            0            0
## 31 2019/2020     2020         366       366            0            0
##    Incomplete_days Perc_complete
## 1                0           100
## 2                0           100
## 3                0           100
## 4                0           100
## 5                0           100
## 6                0           100
## 7                0           100
## 8                0           100
## 9                0           100
## 10               0           100
## 11               0           100
## 12               0           100
## 13               0           100
## 14               0           100
## 15               0           100
## 16               0           100
## 17               0           100
## 18               0           100
## 19               0           100
## 20               0           100
## 21               0           100
## 22               0           100
## 23               0           100
## 24               0           100
## 25               0           100
## 26               0           100
## 27               0           100
## 28               0           100
## 29               0           100
## 30               0           100
## 31               0           100
table(is.na(patched_complete$Tmin))
## 
## FALSE 
## 11323
write.csv(patched_complete, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv", row.names = FALSE)

Chapter 12 Generating temperature scenarios

  1. For the location you chose for your earlier analyses, use chillR’s weather generator to produce 100 years of synthetic temperature data.
patched_complete <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")

# we simulate here with the data form 1998 to 2009 100 times how the weather could have been in this time
Temp <- patched_complete %>%
  temperature_generation(years = c(1998,2009),
                         sim_years = c(2001,2100))
## Warning: scenario doesn't contain named elements - consider using the following
## element names: 'data', 'reference_year','scenario_type','labels'
## Warning: setting scenario_type to 'relative'
## Warning: Reference year missing - can't check if relative temperature scenario
## is valid
Temperatures <- patched_complete %>% filter(Year %in% 1998:2009) %>%
  cbind(Data_source = "observed") # cbind() adds more coloumns to the table

new_data <- Temp[[1]] %>%
  select(Year, Month, Day, Tmin, Tmax) %>%  # Nur die gewünschten Spalten auswählen
  mutate(Data_source = "simulated")        # Neue Spalte Data_source hinzufügen

# Daten an Temperatures anhängen
Temperatures <- bind_rows(Temperatures, new_data)


Temperatures <- Temperatures %>% 
  mutate(Date = as.Date(ISOdate(2000, Month, Day))) # mutate creates a coloumn from already existing data
ggplot(data = Temperatures,
       aes(Date,
           Tmin)) +
  geom_smooth(aes(colour = factor(Year))) +
  facet_wrap(vars(Data_source)) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none") +
  scale_x_date(date_labels = "%b")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(data = Temperatures,
       aes(Date,
           Tmax)) + # creates a empty coordinate system
  geom_smooth(aes(colour = factor(Year))) + # creates a smoothing regression to make the data better obsevable (alternative is geom_line() then you getalso the nice colors but without the smoothing evect)
  facet_wrap(vars(Data_source)) + # 
  theme_bw(base_size = 20) +
  theme(legend.position = "none") + # deletes the legend (is here only confusing)
  scale_x_date(date_labels = "%b") # shows on the xachsis only the month without the year
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

  1. Calculate winter chill (in Chill Portions) for your synthetic weather, and illustrate your results as histograms and cumulative distributions.
chill_observed <- Temperatures %>%
  filter(Data_source == "observed") %>%
  stack_hourly_temps(latitude = 50.4) %>%
  chilling(Start_JDay = 305,
           End_JDay = 59)
  
chill_simulated <- Temperatures %>%
  filter(Data_source == "simulated") %>%
  stack_hourly_temps(latitude = 51.4) %>%
  chilling(Start_JDay = 305,
           End_JDay = 59)
  
chill_comparison <-
  cbind(chill_observed,
        Data_source = "observed") %>%
  rbind(cbind(chill_simulated,
              Data_source = "simulated"))
  
chill_comparison_full_seasons <- 
  chill_comparison %>%
  filter(Perc_complete == 100)

chill_comparison_full_seasons$Data_source <- factor(
  chill_comparison_full_seasons$Data_source,
  levels = c("simulated", "observed")
)


chill_comparison_full_seasons$Data_source <- factor(
  chill_comparison_full_seasons$Data_source,
  levels = c("simulated", "observed")
)

chill_comparison_full_seasons <- chill_comparison_full_seasons[order(
  chill_comparison_full_seasons$Data_source), ]
ggplot(chill_comparison_full_seasons,
       aes(x = Chill_portions)) + 
  geom_histogram(binwidth = 1, position = "identity",
                 aes(fill = factor(Data_source))) +
  scale_fill_brewer(palette="Set2") +
  theme_bw(base_size = 20) +
  labs(fill = "Data source") +
  xlab("Chill accumulation (Chill Portions)") +
  ylab("Frequency")

chill_simulations <-
  chill_comparison_full_seasons %>%
  filter(Data_source == "simulated")
  
ggplot(chill_simulations,
       aes(x = Chill_portions)) +
  stat_ecdf(geom = "step",
            lwd = 1.5,
            col = "blue") +
  ylab("Cumulative probability") +
  xlab("Chill accumulation (in Chill Portions)") +
  theme_bw(base_size = 20)

  1. Produce similar plots for the number of freezing hours (<0°C) in April (or October, if your site is in the Southern Hemisphere) for your location of interest.
# filter only our observed data and calculate with the coordinates (latitude) the chilling day 
chill_observed <- Temperatures %>%
  filter(Data_source == "observed") %>%
  stack_hourly_temps(latitude = 51.4) %>%
  chilling(Start_JDay = 305,
           End_JDay = 59)


### adding a freezing model - to make this for April, you'll also have to adjust the dates in the calculations
df <- data.frame(
  lower= c(-1000, 0),
  upper= c(    0, 1000),
  weight=c(    1, 0))

freezing_hours <- function(x) step_model(x,df)

freezing_hours(c(1, 2, 4, 5, -10))
## [1] 0 0 0 0 1
chill_observed <- Temperatures %>%
  filter(Data_source == "observed") %>%
  stack_hourly_temps(latitude = 51.4) %>%
  tempResponse(Start_JDay = 305, # more flexible than the chilling() function, we feed it our freezing function
               End_JDay = 59,
               models = list(Frost = freezing_hours,
                             Chill_portions = Dynamic_Model,
                             GDH = GDH))


####

# repeating the freezing function with the simulated data
chill_simulated <- Temperatures %>%
  filter(Data_source == "simulated") %>%
  stack_hourly_temps(latitude = 51.4) %>%
  tempResponse(Start_JDay = 305,
               End_JDay = 59,
               models=list(Frost = freezing_hours,
                           Chill_portions = Dynamic_Model,
                           GDH = GDH))

# combine the freezing data for the observed and the simulated
chill_comparison <-
  cbind(chill_observed,
        Data_source = "observed") %>%
  rbind(cbind(chill_simulated,
              Data_source = "simulated"))

# sometimes there is data missing for the winter and the influence the data negative, so we remove them
chill_comparison_full_seasons <-
  chill_comparison %>%
  filter(Perc_complete == 100)


chill_comparison_full_seasons$Data_source <- factor(
  chill_comparison_full_seasons$Data_source,
  levels = c("simulated", "observed")
)

chill_comparison_full_seasons <- chill_comparison_full_seasons[order(
  chill_comparison_full_seasons$Data_source), ]
# plot the data as a histogram
ggplot(chill_comparison_full_seasons,
       aes(x = Frost)) + 
  geom_histogram(binwidth = 25, position = "identity",
                 aes(fill = factor(Data_source))) +
  scale_fill_brewer(palette="Set2") +
  theme_bw(base_size = 10) +
  labs(fill = "Data source") +
  xlab("Frost incidence during winter (hours)") +
  ylab("Frequency")

chill_simulations <-
  chill_comparison_full_seasons %>%
  filter(Data_source == "simulated")

# shows how many frost hours are how likely to happen
ggplot(chill_simulations,
       aes(x = Frost)) +
  stat_ecdf(geom = "step", # accumulated density function
            lwd = 1.5,
            col = "blue") +
  ylab("Cumulative probability") +
  xlab("Frost incidence during winter (hours)") +
  theme_bw(base_size = 20)

# Here's the amount of chill that is exceeded in 90% of all years.
quantile(chill_simulations$Chill_portions, 0.1)
##      10% 
## 78.80063
# and here's the 50% confidence interval (25th to 75th percentile)
quantile(chill_simulations$Chill_portions, c(0.25, 0.75))
##      25%      75% 
## 82.04394 86.97673

Chapter 13 Saving and loading data (and hiding this in markdown)

There are no tasks for this chapter.

Chapter 14 Historic temperature scenarios

  1. For the location you chose for previous exercises, produce historic temperature scenarios representing several years of the historic record (your choice).
patched_complete <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")

Temp <- patched_complete %>%
  temperature_generation(years = c(1998,2005),
                         sim_years = c(2001,2100))

change_scenario <- data.frame(Tmin = rep(2,12),
                              Tmax = rep(2,12))
# change_scenario

Temp_2 <- temperature_generation(KA_weather,
                                 years = c(1998,2005),
                                 sim_years = c(2001,2100),
                                 temperature_scenario = change_scenario)

Temperature_scenarios <- KA_weather %>%
  filter(Year %in% 1998:2005) %>%
  cbind(Data_source = "observed") %>%
  rbind(Temp[[1]] %>% 
          select(c(Year, Month, Day, Tmin, Tmax)) %>% 
          cbind(Data_source = "simulated")
        ) %>%
  rbind(Temp_2[[1]] %>%
          select(c(Year, Month, Day, Tmin, Tmax)) %>% 
          cbind(Data_source = "Warming_2C")
        ) %>%
  mutate(Date = as.Date(ISOdate(2000,
                                Month,
                                Day)))

ggplot(data = Temperature_scenarios, 
       aes(Date,Tmin)) +
  geom_smooth(aes(colour = factor(Year))) +
  facet_wrap(vars(Data_source)) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none") +
  scale_x_date(date_labels = "%b")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

ggplot(data = Temperature_scenarios, 
       aes(Date,Tmax)) +
  geom_smooth(aes(colour = factor(Year))) +
  facet_wrap(vars(Data_source)) +
  theme_bw(base_size = 20) +
  theme(legend.position = "none") +
  scale_x_date(date_labels = "%b")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

summary(patched_complete)
##     YEARMODA            DATE               Date                Year     
##  Min.   :19900101   Length:11323       Length:11323       Min.   :1990  
##  1st Qu.:19971002   Class :character   Class :character   1st Qu.:1997  
##  Median :20050702   Mode  :character   Mode  :character   Median :2005  
##  Mean   :20050675                                         Mean   :2005  
##  3rd Qu.:20130402                                         3rd Qu.:2013  
##  Max.   :20201231                                         Max.   :2020  
##                                                                         
##      Month             Day             Tmin              Tmax        
##  Min.   : 1.000   Min.   : 1.00   Min.   :-19.889   Min.   :-10.278  
##  1st Qu.: 4.000   1st Qu.: 8.00   1st Qu.:  2.276   1st Qu.:  9.222  
##  Median : 7.000   Median :16.00   Median :  6.722   Median : 15.111  
##  Mean   : 6.523   Mean   :15.73   Mean   :  6.484   Mean   : 15.112  
##  3rd Qu.:10.000   3rd Qu.:23.00   3rd Qu.: 11.111   3rd Qu.: 21.000  
##  Max.   :12.000   Max.   :31.00   Max.   : 23.278   Max.   : 40.222  
##                                                                      
##      Tmean              Prec        Tmin_source        Tmax_source       
##  Min.   :-15.222   Min.   : 0.000   Length:11323       Length:11323      
##  1st Qu.:  5.944   1st Qu.: 0.000   Class :character   Class :character  
##  Median : 10.916   Median : 0.000   Mode  :character   Mode  :character  
##  Mean   : 10.784   Mean   : 1.952                                        
##  3rd Qu.: 15.944   3rd Qu.: 2.032                                        
##  Max.   : 31.278   Max.   :78.486                                        
##  NA's   :2301      NA's   :2301                                          
##  Tmin_interpolated Tmax_interpolated
##  Mode :logical     Mode :logical    
##  FALSE:11322       FALSE:11322      
##  TRUE :1           TRUE :1          
##                                     
##                                     
##                                     
## 
scenario_1990 <- temperature_scenario_from_records(weather = patched_complete,
                                                   year = 1990)

temps_1990 <- temperature_generation(weather = patched_complete,
                                     years = c(1990,2020),
                                     sim_years = c(2001,2100),
                                     temperature_scenario = scenario_1990)

scenario_2005 <- temperature_scenario_from_records(weather = patched_complete,
                                                   year = 2005)

relative_scenario <- temperature_scenario_baseline_adjustment(
  baseline = scenario_2005,
  temperature_scenario = scenario_1990)

temps_1990<-temperature_generation(weather = patched_complete,
                                   years = c(1990,2020),
                                   sim_years = c(2001,2100),
                                   temperature_scenario = relative_scenario)

all_past_scenarios <- temperature_scenario_from_records(
  weather = patched_complete,
  year = c(1990,
           1997,
           2003,
           2010))

adjusted_scenarios <- temperature_scenario_baseline_adjustment(
  baseline = scenario_2005,
  temperature_scenario = all_past_scenarios)

all_past_scenario_temps <- temperature_generation(
  weather = patched_complete,
  years = c(1990,2020),
  sim_years = c(2001,2100),
  temperature_scenario = adjusted_scenarios)

save_temperature_scenarios(all_past_scenario_temps, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_scenarios")
  1. Produce chill distributions for these scenarios and plot them.
frost_model <- function(x)
  step_model(x,
             data.frame(
               lower=c(-1000,0),
               upper=c(0,1000),
               weight=c(1,0)))

models <- list(Chill_Portions = Dynamic_Model,
               GDH = GDH,
               Frost_H = frost_model)

chill_hist_scenario_list <- tempResponse_daily_list(all_past_scenario_temps,
                                                    latitude = 51.4,
                                                    Start_JDay = 305,
                                                    End_JDay = 59,
                                                    models = models)

chill_hist_scenario_list <- lapply(chill_hist_scenario_list,
                                   function(x) x %>%
                                     filter(Perc_complete == 100))

save_temperature_scenarios(chill_hist_scenario_list, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_chill_305_59")
scenarios <- names(chill_hist_scenario_list)[1:4]

all_scenarios <- chill_hist_scenario_list[[scenarios[1]]] %>%
  mutate(scenario = as.numeric(scenarios[1]))

for (sc in scenarios[2:4])
 all_scenarios <- all_scenarios %>%
  rbind(chill_hist_scenario_list[[sc]] %>%
          cbind(
            scenario=as.numeric(sc))
        ) %>%
  filter(Perc_complete == 100)

# Let's compute the actual 'observed' chill for comparison
actual_chill <- tempResponse_daily_list(patched_complete,
                                        latitude=51.4,
                                        Start_JDay = 305,
                                        End_JDay = 59,
                                        models)[[1]] %>%
  filter(Perc_complete == 100)

ggplot(data = all_scenarios,
       aes(scenario,
           Chill_Portions,
           fill = factor(scenario))) +
  geom_violin() +
  ylab("Chill accumulation (Chill Portions)") +
  xlab("Scenario year") +
  theme_bw(base_size = 15) +
  ylim(c(0,90)) +
  geom_point(data = actual_chill,
             aes(End_year,
                 Chill_Portions,
                 fill = "blue"),
             col = "blue",
             show.legend = FALSE) +
  scale_fill_discrete(name = "Scenario",
                      breaks = unique(all_scenarios$scenario))
## Warning: Removed 10 rows containing non-finite outside the scale range
## (`stat_ydensity()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

write.csv(actual_chill,"D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_observed_chill_305_59.csv", row.names = FALSE)
temperature_means <- 
  data.frame(Year = min(patched_complete$Year):max(patched_complete$Year),
             Tmin = aggregate(patched_complete$Tmin,
                              FUN = "mean",
                              by = list(patched_complete$Year))[,2],
             Tmax=aggregate(patched_complete$Tmax,
                            FUN = "mean",
                            by = list(patched_complete$Year))[,2]) %>%
  mutate(runn_mean_Tmin = runn_mean(Tmin,15),
         runn_mean_Tmax = runn_mean(Tmax,15))


Tmin_regression <- lm(Tmin~Year,
                      temperature_means)

Tmax_regression <- lm(Tmax~Year,
                      temperature_means)

temperature_means <- temperature_means %>%
  mutate(regression_Tmin = Tmin_regression$coefficients[1]+
           Tmin_regression$coefficients[2]*temperature_means$Year,
           regression_Tmax = Tmax_regression$coefficients[1]+
           Tmax_regression$coefficients[2]*temperature_means$Year
         )

ggplot(temperature_means,
       aes(Year,
           Tmin)) + 
  geom_point() + 
  geom_line(data = temperature_means,
            aes(Year,
                runn_mean_Tmin),
            lwd = 2,
            col = "blue") + 
  geom_line(data = temperature_means,
            aes(Year,
                regression_Tmin),
            lwd = 2,
            col = "red") +
  theme_bw(base_size = 15) +
  ylab("Mean monthly minimum temperature (°C)")

ggplot(temperature_means,
       aes(Year,
           Tmax)) + 
  geom_point() + 
  geom_line(data = temperature_means,
            aes(Year,
                runn_mean_Tmax),
            lwd = 2,
            col = "blue") + 
  geom_line(data = temperature_means,
            aes(Year, 
                regression_Tmax),
            lwd = 2,
            col = "red") +
  theme_bw(base_size = 15) +
  ylab("Mean monthly maximum temperature (°C)")

Chapter 15 Future temperature scenarios

  1. Briefly describe the differences between the RCPs and the SSPs (you may have to follow some of the links provided above).

The main differences between RCPs (Representative Concentration Pathways) ans SSPs (Shared Socioeconomic Pathways) are the following:

  • Focus

    • RCPs are defined through greenhouse gas concentrations and radioactive forcing levels

    • SSPs include socioeconomic factors (e.g. population growth, energy use, policy decisions) along with emissions pathways

  • Drivers of Change

    • RCPs only consider atmospheric greenhouse gas concentrations and their effects on climate

    • SSPs include also social, economic and technological factors that influence emissions

  • Policy Relevance

    • RCPs provide no direct link to human actions or policy decisions

    • SSPs offer insights into how different societal choices impact climate change

  • Emissions Pathway Representation

    • RCPs directly define future greenhouse gas concentration levels

    • SSPs can be paired with radioactive forcing levels for more flexibility

  • RCPs are considered now as outdated and SSPs are now mainly used

Chapter 16 Making CMIP6 scenarios

  1. Analyse the historic and future impact of climate change on two agroclimatic metrics of your choice, for the location you’ve chosen for your earlier analyses.
location = c(6.29, 51.40)
area <- c(52.5, 5, 50.5, 7)

download_cmip6_ecmwfr(
  scenarios = c("ssp126", "ssp245", "ssp370", "ssp585"),
  area = area,
  key = '9d522e37-4b90-46d1-8701-7deadd45032b',
  model = 'default',
  frequency = 'monthly',
  variable = c('Tmin', 'Tmax'),
  year_start = 2015,
  year_end = 2100)
download_baseline_cmip6_ecmwfr(
  area = area,
  key = '9d522e37-4b90-46d1-8701-7deadd45032b',
  model = 'match_downloaded',
  frequency = 'monthly',
  variable = c('Tmin', 'Tmax'),
  year_start = 1986,
  year_end = 2014,
  month = 1:12)
library(ncdf4)
library(PCICt)
Wesel_temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")

station <- data.frame(
  station_name = c("ARCEN AWS"),
  longitude = c(6.5),
  latitude = c(51.8))

extracted <- extract_cmip6_data(stations = station)
## Unzipping files
## Extracting downloaded CMIP6 files
# head(extracted$`ssp126_AWI-CM-1-1-MR`)

change_scenarios <- 
  gen_rel_change_scenario(extracted,
                          scenarios = c(2050, 2085),
                          reference_period = c(1986:2014),
                          future_window_width = 30)
head(change_scenarios)
## # A tibble: 6 × 11
##   location  Month  Tmax  Tmin scenario start_year end_year scenario_year
##   <chr>     <dbl> <dbl> <dbl> <chr>         <dbl>    <dbl>         <dbl>
## 1 ARCEN AWS     1  1.90  1.95 ssp126         2035     2065          2050
## 2 ARCEN AWS     2  2.60  2.00 ssp126         2035     2065          2050
## 3 ARCEN AWS     3  1.40  1.04 ssp126         2035     2065          2050
## 4 ARCEN AWS     4  1.38  1.15 ssp126         2035     2065          2050
## 5 ARCEN AWS     5  1.87  1.56 ssp126         2035     2065          2050
## 6 ARCEN AWS     6  1.86  1.67 ssp126         2035     2065          2050
## # ℹ 3 more variables: reference_year <int>, scenario_type <chr>, labels <chr>
write.csv(change_scenarios, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/all_change_scenarios_ssp126.csv", row.names = FALSE)
change_scenarios <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/all_change_scenarios_ssp126.csv")


scen_list <- convert_scen_information(change_scenarios)

scen_frame <- convert_scen_information(scen_list)

# scen_list$"ARCEN AWS"$ssp126$`ACCESS-CM2`$'2050'
temps_2005 <- temperature_scenario_from_records(Wesel_temps,
                                                2005)
temps_2000 <- temperature_scenario_from_records(Wesel_temps,
                                                2000)
# temps_2005
# temps_2000

# Wesel_temps
base <- temperature_scenario_baseline_adjustment(temps_2005,
                                                 temps_2000)

scen_list <- convert_scen_information(change_scenarios,
                                      give_structure = FALSE)

adjusted_list <- temperature_scenario_baseline_adjustment(base,
                                                          scen_list,
                   temperature_check_args = 
                     list(scenario_check_thresholds = c(-5, 15)))
temps <- temperature_generation(Wesel_temps,
                                years = c(1990, 2020),
                                sim_years = c(2001, 2100),
                                adjusted_list,  
                                temperature_check_args = 
                                  list( scenario_check_thresholds = c(-5, 15)))

save_temperature_scenarios(temps, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_future")
for(scen in 1:length(adjusted_list))
{
  if(!file.exists(paste0("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate/Wesel_future_",
                       scen,"_",
                       names(adjusted_list)[scen],".csv")) )
  {temp_temp <- temperature_generation(Wesel_temps,
                                   years = c(1990, 2020),
                                   sim_years = c(2001, 2100),
                                   adjusted_list[scen],  
                                   temperature_check_args = 
                                     list( scenario_check_thresholds = c(-5, 15)))
  write.csv(temp_temp[[1]],paste0("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate/Wesel_future_",scen,"_",names(adjusted_list)[scen],".csv"),
                             row.names=FALSE)
  print(paste("Processed object",scen,"of", length(adjusted_list)))
  
    
  }
  
}
temps <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_future_")

frost_model <- function(x)
  step_model(x,
             data.frame(
               lower = c(-1000, 0),
               upper = c(0, 1000),
               weight = c(1, 0)))

models <- list(Chill_Portions = Dynamic_Model,
               GDH = GDH,
               Frost_H = frost_model)
chill_future_scenario_list <- 
  tempResponse_daily_list(temps,
                          latitude = 51.4,
                          Start_JDay = 305,
                          End_JDay = 59,
                          models = models)

chill_future_scenario_list <- 
  lapply(chill_future_scenario_list,
         function(x) x %>%
           filter(Perc_complete == 100))

save_temperature_scenarios(chill_future_scenario_list,
                           "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
                           "Wesel_futurechill")

chill_future_scenario_list_frost <- 
  lapply(chill_future_scenario_list_frost,
         function(x) x %>%
           filter(Perc_complete == 100))


save_temperature_scenarios(chill_future_scenario_list_frost,
                           "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
                           "Wesel_futurefrost")
Wesel_temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")

observed_chill <- tempResponse_daily_list(Wesel_temps,
                                          latitude = 51.4,
                                          Start_JDay = 305,
                                          End_JDay = 59,
                                          models = models)
observed_chill <- lapply(observed_chill,
                         function(x) x %>%
                           filter(Perc_complete == 100))

write.csv(observed_chill, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_observed_chill_305_59.csv")
hist_temps <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_chill_")

chill_hist_scenario_list <- tempResponse_daily_list(hist_temps,
                                          latitude = 51.4,
                                          Start_JDay = 305,
                                          End_JDay = 59,
                                          models = models)

chill_hist_scenario_list <- lapply(chill_hist_scenario_list,
                                     function(x) x %>%
                                       filter(Perc_complete == 100))


save_temperature_scenarios(chill_hist_scenario_list,
                           "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
                           "Wesel_hist_chill_305_59")
chill_future_scenario_list <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_futurechill")

chill_hist_scenario_list<-load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_chill_305_59")

observed_chill <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_observed_chill_305_59.csv")

chills <- make_climate_scenario(
  metric_summary = chill_hist_scenario_list,
  caption = "Historical",
  historic_data = observed_chill,
  time_series = TRUE)

plot_climate_scenarios(
  climate_scenario_list = chills,
  metric = "Chill_Portions",
  metric_label = "Chill (Chill Portions)")

## [[1]]
## [1] "time series labels"
list_ssp <- 
  strsplit(names(chill_future_scenario_list),
           '\\.') %>%
  map(2) %>%
  unlist()

list_gcm <-
  strsplit(names(chill_future_scenario_list), '\\.') %>%
  map(3) %>%
  unlist()

list_time <-
  strsplit(names(chill_future_scenario_list), '\\.') %>%
  map(4) %>%
  unlist()

SSPs <- c("ssp126", "ssp245", "ssp370", "ssp585")
Times <- c(2050, 2085)


for(SSP in SSPs)
  for(Time in Times)
    {
    
    # find all scenarios for the ssp and time
    chill <- chill_future_scenario_list[list_ssp == SSP &
                                          list_time == Time]
    names(chill) <- list_gcm[list_ssp == SSP &
                               list_time == Time]
    if(SSP == "ssp126") SSPcaption <- "SSP1"
    if(SSP == "ssp245") SSPcaption <- "SSP2"
    if(SSP == "ssp370") SSPcaption <- "SSP3"
    if(SSP == "ssp585") SSPcaption <- "SSP5"    
    chills <- chill %>% 
      make_climate_scenario(
        caption = c(SSPcaption,
                    Time),
        add_to = chills)
}

info_chill <-
  plot_climate_scenarios(
    climate_scenario_list = chills,
    metric = "Chill_Portions",
    metric_label = "Chill (Chill Portions)",
    texcex = 1)

info_heat <-
  plot_climate_scenarios(
    climate_scenario_list = chills,
    metric = "GDH",
    metric_label = "Heat (Growing Degree Hours)",
    texcex = 1)

info_frost <- 
  plot_climate_scenarios(  
    climate_scenario_list=chills,
    metric="Frost_H",
    metric_label="Frost incidence (hours)",
    texcex=1)

info_chill[[2]]
##    code            Label
## 1     1        INM-CM4-8
## 2     2    MPI-ESM1-2-LR
## 3     3        CMCC-ESM2
## 4     4       MIROC-ES2L
## 5     5      FIO-ESM-2-0
## 6     6    CNRM-CM6-1-HR
## 7     7        INM-CM5-0
## 8     8       MRI-ESM2-0
## 9     9 EC-Earth3-Veg-LR
## 10   10        GFDL-ESM4
## 11   11           MIROC6
## 12   12      CNRM-ESM2-1
## 13   13     IPSL-CM6A-LR
## 14   14            NESM3
## 15   15        FGOALS-g3
## 16   16       ACCESS-CM2
## 17   17    AWI-CM-1-1-MR
## 18   18          CanESM5

Chapter 17 Making CMIP5 scenarios with the ClimateWizard

There are no tasks for this chapter.

Chapter 18 Plotting future scenarios

  1. Produce similar plots for the weather station you selected for earlier exercises.
chill_hist_scenario_list <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_chill_305_59")
actual_chill <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_observed_chill_305_59.csv")

chill_future_scenario_list <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_futurechill")

chills <- make_climate_scenario(
  chill_hist_scenario_list,
  caption = "Historic",
  historic_data = actual_chill,
  time_series = TRUE)

SSPs <- c("ssp126", "ssp245", "ssp370", "ssp585")
Times <- c(2050, 2085)

list_ssp <- 
  strsplit(names(chill_future_scenario_list), '\\.') %>%
  map(2) %>%
  unlist()

list_gcm <-
  strsplit(names(chill_future_scenario_list), '\\.') %>%
  map(3) %>%
  unlist()

list_time <-
  strsplit(names(chill_future_scenario_list), '\\.') %>%
  map(4) %>%
  unlist()


for(SSP in SSPs)
  for(Time in Times)
    {
    
    # find all scenarios for the ssp and time
    chill <- chill_future_scenario_list[list_ssp == SSP & list_time == Time]
    names(chill) <- list_gcm[list_ssp == SSP & list_time == Time]
    if(SSP == "ssp126") SSPcaption <- "SSP1"
    if(SSP == "ssp245") SSPcaption <- "SSP2"
    if(SSP == "ssp370") SSPcaption <- "SSP3"
    if(SSP == "ssp585") SSPcaption <- "SSP5"    
    if(Time == "2050") Time_caption <- "2050"
    if(Time == "2085") Time_caption <- "2085"
    chills <- chill %>% 
      make_climate_scenario(
        caption = c(SSPcaption,
                    Time_caption),
        add_to = chills)
  }

plot_climate_scenarios(
  climate_scenario_list = chills,
  metric = "Chill_Portions",
  metric_label = "Chill (Chill Portions)",
  texcex = 1)

## [[1]]
## [1] "time series labels"
## 
## [[2]]
##    code            Label
## 1     1        INM-CM4-8
## 2     2    MPI-ESM1-2-LR
## 3     3        CMCC-ESM2
## 4     4       MIROC-ES2L
## 5     5      FIO-ESM-2-0
## 6     6    CNRM-CM6-1-HR
## 7     7        INM-CM5-0
## 8     8       MRI-ESM2-0
## 9     9 EC-Earth3-Veg-LR
## 10   10        GFDL-ESM4
## 11   11           MIROC6
## 12   12      CNRM-ESM2-1
## 13   13     IPSL-CM6A-LR
## 14   14            NESM3
## 15   15        FGOALS-g3
## 16   16       ACCESS-CM2
## 17   17    AWI-CM-1-1-MR
## 18   18          CanESM5
## 
## [[3]]
##    code            Label
## 1     1       ACCESS-CM2
## 2     2        CMCC-ESM2
## 3     3      FIO-ESM-2-0
## 4     4    CNRM-CM6-1-HR
## 5     5    AWI-CM-1-1-MR
## 6     6 EC-Earth3-Veg-LR
## 7     7        GFDL-ESM4
## 8     8      CNRM-ESM2-1
## 9     9          CanESM5
## 10   10        FGOALS-g3
## 11   11        INM-CM4-8
## 12   12        INM-CM5-0
## 13   13     IPSL-CM6A-LR
## 14   14       MIROC-ES2L
## 15   15           MIROC6
## 16   16    MPI-ESM1-2-LR
## 17   17       MRI-ESM2-0
## 18   18            NESM3
## 
## [[4]]
##    code            Label
## 1     1        FGOALS-g3
## 2     2        INM-CM4-8
## 3     3    MPI-ESM1-2-LR
## 4     4        CMCC-ESM2
## 5     5     EC-Earth3-CC
## 6     6       MIROC-ES2L
## 7     7      FIO-ESM-2-0
## 8     8    CNRM-CM6-1-HR
## 9     9        INM-CM5-0
## 10   10       MRI-ESM2-0
## 11   11 EC-Earth3-Veg-LR
## 12   12        GFDL-ESM4
## 13   13           MIROC6
## 14   14      CNRM-ESM2-1
## 15   15     IPSL-CM6A-LR
## 16   16            NESM3
## 17   17       ACCESS-CM2
## 18   18    AWI-CM-1-1-MR
## 
## [[5]]
##    code            Label
## 1     1        FGOALS-g3
## 2     2       ACCESS-CM2
## 3     3        CMCC-ESM2
## 4     4     EC-Earth3-CC
## 5     5      FIO-ESM-2-0
## 6     6    CNRM-CM6-1-HR
## 7     7    AWI-CM-1-1-MR
## 8     8 EC-Earth3-Veg-LR
## 9     9        GFDL-ESM4
## 10   10      CNRM-ESM2-1
## 11   11        INM-CM4-8
## 12   12        INM-CM5-0
## 13   13     IPSL-CM6A-LR
## 14   14       MIROC-ES2L
## 15   15           MIROC6
## 16   16    MPI-ESM1-2-LR
## 17   17       MRI-ESM2-0
## 18   18            NESM3
## 
## [[6]]
##    code             Label
## 1     1       CNRM-ESM2-1
## 2     2      IPSL-CM6A-LR
## 3     3         FGOALS-g3
## 4     4 EC-Earth3-AerChem
## 5     5         INM-CM4-8
## 6     6     MPI-ESM1-2-LR
## 7     7        MIROC-ES2L
## 8     8     CNRM-CM6-1-HR
## 9     9         INM-CM5-0
## 10   10        MRI-ESM2-0
## 11   11  EC-Earth3-Veg-LR
## 12   12         GFDL-ESM4
## 13   13            MIROC6
## 14   14        ACCESS-CM2
## 15   15     AWI-CM-1-1-MR
## 
## [[7]]
##    code             Label
## 1     1       CNRM-ESM2-1
## 2     2         FGOALS-g3
## 3     3 EC-Earth3-AerChem
## 4     4        ACCESS-CM2
## 5     5     CNRM-CM6-1-HR
## 6     6     AWI-CM-1-1-MR
## 7     7  EC-Earth3-Veg-LR
## 8     8         GFDL-ESM4
## 9     9         INM-CM4-8
## 10   10         INM-CM5-0
## 11   11      IPSL-CM6A-LR
## 12   12        MIROC-ES2L
## 13   13            MIROC6
## 14   14     MPI-ESM1-2-LR
## 15   15        MRI-ESM2-0
## 
## [[8]]
##    code            Label
## 1     1            CIESM
## 2     2            NESM3
## 3     3      CNRM-ESM2-1
## 4     4     IPSL-CM6A-LR
## 5     5        FGOALS-g3
## 6     6        CMCC-ESM2
## 7     7        INM-CM4-8
## 8     8    MPI-ESM1-2-LR
## 9     9     EC-Earth3-CC
## 10   10      FIO-ESM-2-0
## 11   11       MIROC-ES2L
## 12   12    CNRM-CM6-1-HR
## 13   13        INM-CM5-0
## 14   14       MRI-ESM2-0
## 15   15 EC-Earth3-Veg-LR
## 16   16        GFDL-ESM4
## 17   17           MIROC6
## 18   18       ACCESS-CM2
## 19   19    AWI-CM-1-1-MR
## 
## [[9]]
##    code            Label
## 1     1            CIESM
## 2     2      CNRM-ESM2-1
## 3     3        FGOALS-g3
## 4     4        CMCC-ESM2
## 5     5       ACCESS-CM2
## 6     6     EC-Earth3-CC
## 7     7      FIO-ESM-2-0
## 8     8    CNRM-CM6-1-HR
## 9     9    AWI-CM-1-1-MR
## 10   10 EC-Earth3-Veg-LR
## 11   11        GFDL-ESM4
## 12   12        INM-CM4-8
## 13   13        INM-CM5-0
## 14   14     IPSL-CM6A-LR
## 15   15       MIROC-ES2L
## 16   16           MIROC6
## 17   17    MPI-ESM1-2-LR
## 18   18       MRI-ESM2-0
## 19   19            NESM3
for(nam in names(chills[[1]]$data))
  {
   ch <- chills[[1]]$data[[nam]]
   ch[,"GCM"] <- "none"
   ch[,"SSP"] <- "none"
   ch[,"Year"] <- as.numeric(nam)
   
  if(nam == names(chills[[1]]$data)[1])
    past_simulated <- ch else
      past_simulated <- rbind(past_simulated,
                              ch)
  }

past_simulated["Scenario"] <- "Historical"

head(past_simulated)
##      Season End_year Season_days Data_days Perc_complete Chill_Portions
## 1 2001/2002     2002         120       120           100       85.93502
## 2 2002/2003     2003         120       120           100       84.11282
## 3 2003/2004     2004         120       120           100       83.33690
## 4 2004/2005     2005         121       121           100       88.70952
## 5 2005/2006     2006         120       120           100       89.60029
## 6 2006/2007     2007         120       120           100       84.94518
##        GDH Frost_H  GCM  SSP Year   Scenario
## 1 3996.664     450 none none 1990 Historical
## 2 3344.907     479 none none 1990 Historical
## 3 1889.238     721 none none 1990 Historical
## 4 2464.805     361 none none 1990 Historical
## 5 3750.400     214 none none 1990 Historical
## 6 3091.575     432 none none 1990 Historical
past_simulated <- past_simulated %>% filter(Perc_complete == 100)

past_observed <- chills[[1]][["historic_data"]]

head(past_observed)
##   X    Season End_year Season_days Data_days Perc_complete Chill_Portions
## 1 1 1990/1991     1991         120       120           100       77.46451
## 2 2 1991/1992     1992         120       120           100       84.12147
## 3 3 1992/1993     1993         121       121           100       85.89337
## 4 4 1993/1994     1994         120       120           100       83.83349
## 5 5 1994/1995     1995         120       120           100       86.00229
## 6 6 1995/1996     1996         120       120           100       67.05746
##        GDH Frost_H
## 1 1964.757     839
## 2 2481.648     491
## 3 3143.898     482
## 4 1904.884     591
## 5 5425.112     317
## 6 1612.944    1224
for(i in 2:length(chills))
  for(nam in names(chills[[i]]$data))
    {ch <- chills[[i]]$data[[nam]]
     ch[,"GCM"] <- nam
     ch[,"SSP"] <- chills[[i]]$caption[1]
     ch[,"Year"] <- chills[[i]]$caption[2]
     if(i == 2 & nam == names(chills[[i]]$data)[1])
       future_data <- ch else
         future_data <- rbind(future_data,ch)
  }

head(future_data)
##      Season End_year Season_days Data_days Perc_complete Chill_Portions
## 1 2001/2002     2002         120       120           100       86.99319
## 2 2002/2003     2003         120       120           100       86.87197
## 3 2003/2004     2004         120       120           100       85.62412
## 4 2004/2005     2005         121       121           100       90.50165
## 5 2005/2006     2006         120       120           100       90.64776
## 6 2006/2007     2007         120       120           100       86.32296
##        GDH Frost_H       GCM  SSP Year
## 1 5318.501     352 INM-CM4-8 SSP1 2050
## 2 4400.557     409 INM-CM4-8 SSP1 2050
## 3 2744.363     578 INM-CM4-8 SSP1 2050
## 4 3436.639     242 INM-CM4-8 SSP1 2050
## 5 5051.034     111 INM-CM4-8 SSP1 2050
## 6 4095.937     350 INM-CM4-8 SSP1 2050
metric <- "GDH"
axis_label <- "Heat (in GDH)"

rng <- range(past_observed[[metric]],
             past_simulated[[metric]],
             future_data[[metric]])  
rng
## [1]  1088.303 19014.867
past_plot <- ggplot() +
  geom_boxplot(data = past_simulated,
               aes_string("as.numeric(Year)",
                          metric,
                          group = "Year"),
               fill = "skyblue")
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
past_plot <- past_plot +
  scale_y_continuous(
    limits = c(0, 
               round(rng[2] + rng[2]/10))) +
  labs(x = "Year", 
       y = axis_label)

past_plot <- past_plot +
  facet_grid(~ Scenario) +
  theme_bw(base_size = 15) 

past_plot <- past_plot +  
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        axis.text.x = element_text(angle = 45,
                                   hjust = 1)) 

# add historic data
past_plot <- past_plot +
  geom_point(data = past_observed,
             aes_string("End_year",
                        metric),
             col = "blue")

past_plot

y <- 2050

future_2050 <-
  ggplot(data = future_data[future_data$Year == y,]) +
  geom_boxplot(aes_string("GCM", 
                          metric, 
                          fill = "GCM"))

future_2050 <- future_2050 +
  facet_wrap(vars(SSP), nrow = 1) +
   scale_x_discrete(labels = NULL,
                    expand = expansion(add = 1)) 

library(ggpmisc)

future_2050 <- future_2050 +
  scale_y_continuous(limits = c(0, 
                                round(round(1.1*rng[2])))) +
    geom_text_npc(aes(npcx = "center", 
                      npcy = "top",
                      label = Year),
                  size = 5)

future_2050 <- future_2050 +
  theme_bw(base_size = 15) +
  theme(axis.ticks.y = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        legend.position = "bottom",
        legend.margin = margin(0,
                               0,
                               0,
                               0,
                               "cm"),
        legend.background = element_rect(),
        strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        legend.box.spacing = unit(0, "cm"),
        plot.subtitle = element_text(hjust = 0.5,
                                     vjust = -1,
                                     size = 15 * 1.05,
                                     face = "bold")) 

future_2050

future_plot_list <- list()

time_points <- c(2050, 2085)

for(y in time_points)
{
  future_plot_list[[which(y == time_points)]] <-
    ggplot(data = future_data[which(future_data$Year==y),]) +
    geom_boxplot(aes_string("GCM",
                            metric,
                            fill="GCM")) +
    facet_wrap(vars(SSP), nrow = 1) +
    scale_x_discrete(labels = NULL,
                     expand = expansion(add = 1)) +
    scale_y_continuous(limits = c(0, 
                                  round(round(1.1*rng[2])))) +
    geom_text_npc(aes(npcx = "center",
                      npcy = "top", 
                      label = Year),
                  size = 5) +
    theme_bw(base_size = 15) +
    theme(axis.ticks.y = element_blank(),
          axis.text = element_blank(),
          axis.title = element_blank(),
          legend.position = "bottom",
          legend.margin = margin(0, 
                                 0, 
                                 0, 
                                 0, 
                                 "cm"),
          legend.background = element_rect(),
          strip.background = element_blank(),
          strip.text = element_text(face = "bold"),
          legend.box.spacing = unit(0, "cm"),
          plot.subtitle = element_text(
            hjust = 0.5,
            vjust = -1,
            size = 15 * 1.05,
            face = "bold")) 
}

library(patchwork)

both_plots <- past_plot + future_plot_list

plot <- both_plots +
           plot_layout(guides = "collect",
                       widths = c(1,
                                  rep(1.8,
                                      length(future_plot_list))))

plot <- plot & theme(legend.position = "bottom",
                     legend.text = element_text(size = 8),
                     legend.title = element_text(size = 10),
                     axis.title.x = element_blank())

plot

metric <- "Chill_Portions"
axis_label <- "Chill (in CP)"

rng <- range(past_observed[[metric]],
             past_simulated[[metric]],
             future_data[[metric]])  
past_plot <- ggplot() +
  geom_boxplot(data = past_simulated,
               aes_string("as.numeric(Year)",
                          metric,group="Year"),
               fill="skyblue") +
  scale_y_continuous(limits = c(0, 
                                round(round(1.1*rng[2])))) +
  labs(x = "Year", y = axis_label) +
  facet_grid(~ Scenario) +
  theme_bw(base_size = 15) +  
  theme(strip.background = element_blank(),
           strip.text = element_text(face = "bold"),
           axis.text.x = element_text(angle=45, hjust=1)) +
  geom_point(data = past_observed,
             aes_string("End_year",
                        metric),
             col="blue")

future_plot_list <- list()

for(y in c(2050,
           2085))
{
  future_plot_list[[which(y == c(2050,2085))]] <-
    ggplot(data = future_data[which(future_data$Year==y),]) +
    geom_boxplot(aes_string("GCM", 
                            metric, 
                            fill="GCM")) +
  facet_wrap(vars(SSP), nrow = 1) +
   scale_x_discrete(labels = NULL,
                    expand = expansion(add = 1)) +
  scale_y_continuous(limits = c(0,
                                round(round(1.1*rng[2])))) +
    geom_text_npc(aes(npcx = "center", 
                      npcy = "top", 
                      label = Year),
                  size = 5) +
  theme_bw(base_size = 15) +
  theme(axis.ticks.y = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        legend.position = "bottom",
        legend.margin = margin(0,
                               0, 
                               0,
                               0, 
                               "cm"),
        legend.background = element_rect(),
        strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        legend.box.spacing = unit(0, "cm"),
        plot.subtitle = element_text(hjust = 0.5, 
                                     vjust = -1,
                                     size = 15 * 1.05,
                                     face = "bold")) 
}

plot <- (past_plot +
           future_plot_list +
           plot_layout(guides = "collect",
                       widths = c(1,rep(1.8,length(future_plot_list))))
         ) & theme(legend.position = "bottom",
                   legend.text = element_text(size = 8),
                   legend.title = element_text(size = 10),
                   axis.title.x=element_blank())
plot

metric <- "Frost_H"
axis_label <- "Frost duration (in hours)"

# get extreme values for the axis scale

rng <- range(past_observed[[metric]],
             past_simulated[[metric]],
             future_data[[metric]])  
past_plot <- ggplot() +
  geom_boxplot(data = past_simulated,
               aes_string("as.numeric(Year)",
                          metric,group="Year"),
               fill="skyblue") +
  scale_y_continuous(limits = c(0, 
                                round(round(1.1*rng[2])))) +
  labs(x = "Year", y = axis_label) +
  facet_grid(~ Scenario) +
  theme_bw(base_size = 15) +  
  theme(strip.background = element_blank(),
           strip.text = element_text(face = "bold"),
           axis.text.x = element_text(angle=45, hjust=1)) +
  geom_point(data = past_observed,
             aes_string("End_year",
                        metric),
             col="blue")

future_plot_list <- list()

for(y in c(2050,
           2085))
{
  future_plot_list[[which(y == c(2050,2085))]] <-
    ggplot(data = future_data[which(future_data$Year==y),]) +
    geom_boxplot(aes_string("GCM", 
                            metric, 
                            fill="GCM")) +
  facet_wrap(vars(SSP), nrow = 1) +
   scale_x_discrete(labels = NULL,
                    expand = expansion(add = 1)) +
  scale_y_continuous(limits = c(0,
                                round(round(1.1*rng[2])))) +
    geom_text_npc(aes(npcx = "center", 
                      npcy = "top", 
                      label = Year),
                  size = 5) +
  theme_bw(base_size = 15) +
  theme(axis.ticks.y = element_blank(),
        axis.text = element_blank(),
        axis.title = element_blank(),
        legend.position = "bottom",
        legend.margin = margin(0,
                               0, 
                               0,
                               0, 
                               "cm"),
        legend.background = element_rect(),
        strip.background = element_blank(),
        strip.text = element_text(face = "bold"),
        legend.box.spacing = unit(0, "cm"),
        plot.subtitle = element_text(hjust = 0.5, 
                                     vjust = -1,
                                     size = 15 * 1.05,
                                     face = "bold")) 
}

plot <- (past_plot +
           future_plot_list +
           plot_layout(guides = "collect",
                       widths = c(1,rep(1.8,length(future_plot_list))))
         ) & theme(legend.position = "bottom",
                   legend.text = element_text(size = 8),
                   legend.title = element_text(size = 10),
                   axis.title.x=element_blank())
plot

plot_scenarios_gg <- function(past_observed,
                              past_simulated,
                              future_data,
                              metric,
                              axis_label,
                              time_points)
{
  rng <- range(past_observed[[metric]],
               past_simulated[[metric]],
               future_data[[metric]])  
  past_plot <- ggplot() +
    geom_boxplot(data = past_simulated,
                 aes_string("as.numeric(Year)",
                            metric,
                            group="Year"),
                 fill="skyblue") +
    scale_y_continuous(limits = c(0, 
                                  round(round(1.1*rng[2])))) +
    labs(x = "Year", y = axis_label) +
    facet_grid(~ Scenario) +
    theme_bw(base_size = 15) +  
    theme(strip.background = element_blank(),
          strip.text = element_text(face = "bold"),
          axis.text.x = element_text(angle=45, 
                                     hjust=1)) +
    geom_point(data = past_observed,
               aes_string("End_year",
                          metric),
               col="blue")
  
  future_plot_list <- list()
  
  for(y in time_points)
  {
    future_plot_list[[which(y == time_points)]] <-
      ggplot(data = future_data[which(future_data$Year==y),]) +
      geom_boxplot(aes_string("GCM", 
                              metric, 
                              fill="GCM")) +
      facet_wrap(vars(SSP), nrow = 1) +
      scale_x_discrete(labels = NULL,
                       expand = expansion(add = 1)) +
      scale_y_continuous(limits = c(0, 
                                    round(round(1.1*rng[2])))) +
      geom_text_npc(aes(npcx = "center",
                        npcy = "top",
                        label = Year),
                    size = 5) +
      theme_bw(base_size = 15) +
      theme(axis.ticks.y = element_blank(),
            axis.text = element_blank(),
            axis.title = element_blank(),
            legend.position = "bottom",
            legend.margin = margin(0,
                                   0, 
                                   0, 
                                   0, 
                                   "cm"),
            legend.background = element_rect(),
            strip.background = element_blank(),
            strip.text = element_text(face = "bold"),
            legend.box.spacing = unit(0, "cm"),
            plot.subtitle = element_text(hjust = 0.5,
                                         vjust = -1,
                                         size = 15 * 1.05,
                                         face = "bold")) 
  }
  
  plot <- (past_plot +
             future_plot_list +
             plot_layout(guides = "collect",
                         widths = c(1,rep(1.8,length(future_plot_list))))
           ) & theme(legend.position = "bottom",
                     legend.text = element_text(size=8),
                     legend.title = element_text(size=10),
                     axis.title.x=element_blank())
  plot
  
}
plot_scenarios_gg(past_observed = past_observed,
                  past_simulated = past_simulated,
                  future_data = future_data,
                  metric = "GDH",
                  axis_label = "Heat (in Growing Degree Hours)",
                  time_points = c(2050, 2085))

ggsave("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_heat_plot.png", 
       width = 10,
       height = 8,
       dpi = 600)

plot_scenarios_gg(past_observed = past_observed,
                  past_simulated = past_simulated,
                  future_data = future_data,
                  metric = "Chill_Portions",
                  axis_label = "Chill (in Chill Portions)",
                  time_points = c(2050, 2085))

ggsave("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_Chill_portions_plot.png", 
       width = 10,
       height = 8,
       dpi = 600)

plot_scenarios_gg(past_observed = past_observed,
                  past_simulated = past_simulated,
                  future_data = future_data,
                  metric = "Frost_H",
                  axis_label = "Frost duration (in hours)",
                  time_points = c(2050, 2085))

ggsave("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Wesel_frost_plot.png", 
       width = 10,
       height = 8,
       dpi = 600)

Chapter 19 Chill model comparison

  1. Perform a similar analysis for the location you’ve chosen for your exercises.
hourly_models <- list(Chilling_units = chilling_units,
     Low_chill = low_chill_model,
     Modified_Utah = modified_utah_model,
     North_Carolina = north_carolina_model,
     Positive_Utah = positive_utah_model,
     Chilling_Hours = Chilling_Hours,
     Utah_Chill_Units = Utah_Model,
     Chill_Portions = Dynamic_Model)

daily_models <- list(Rate_of_Chill = rate_of_chill,
                     Chill_Days = chill_days,
                     Exponential_Chill = exponential_chill,
                     # Triangular_Chill_Haninnen = triangular_chill_1,
                     Triangular_Chill_Legave = triangular_chill_2)

metrics <- c(names(daily_models),
             names(hourly_models))

model_labels = c("Rate of Chill",
                 "Chill Days",
                 "Exponential Chill",
                 # "Triangular Chill (Häninnen)",
                 "Triangular Chill (Legave)",
                 "Chilling Units",
                 "Low-Chill Chill Units",
                 "Modified Utah Chill Units",
                 "North Carolina Chill Units",
                 "Positive Utah Chill Units",
                 "Chilling Hours",
                 "Utah Chill Units",
                 "Chill Portions")

data.frame(Metric = model_labels, 'Function name' = metrics)
##                        Metric           Function.name
## 1               Rate of Chill           Rate_of_Chill
## 2                  Chill Days              Chill_Days
## 3           Exponential Chill       Exponential_Chill
## 4   Triangular Chill (Legave) Triangular_Chill_Legave
## 5              Chilling Units          Chilling_units
## 6       Low-Chill Chill Units               Low_chill
## 7   Modified Utah Chill Units           Modified_Utah
## 8  North Carolina Chill Units          North_Carolina
## 9   Positive Utah Chill Units           Positive_Utah
## 10             Chilling Hours          Chilling_Hours
## 11           Utah Chill Units        Utah_Chill_Units
## 12             Chill Portions          Chill_Portions
Wesel_temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv")

Temps <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters", "Wesel_hist_scenarios")
Start_JDay <- 305
End_JDay <- 59

daily_models_past_scenarios <- 
  tempResponse_list_daily(Temps,
                          Start_JDay = Start_JDay,
                          End_JDay = End_JDay,
                          models=daily_models)

daily_models_past_scenarios <- lapply(
  daily_models_past_scenarios,
  function(x) x[which(x$Perc_complete>90),])

hourly_models_past_scenarios<-
  tempResponse_daily_list(Temps,
                          latitude = 51.405,
                          Start_JDay = Start_JDay,
                          End_JDay = End_JDay,
                          models = hourly_models,
                          misstolerance = 10)

past_scenarios <- daily_models_past_scenarios

past_scenarios <- lapply(
  names(past_scenarios),
  function(x)
    cbind(past_scenarios[[x]],
          hourly_models_past_scenarios[[x]][,names(hourly_models)]))

names(past_scenarios) <- names(daily_models_past_scenarios)

daily_models_observed <- 
  tempResponse_daily(Wesel_temps,
                     Start_JDay = Start_JDay,
                     End_JDay = End_JDay,
                     models = daily_models)

daily_models_observed <-
  daily_models_observed[which(daily_models_observed$Perc_complete>90),]

hourly_models_observed <- 
  tempResponse_daily_list(Wesel_temps,
                          latitude=51.405,
                          Start_JDay = Start_JDay,
                          End_JDay = End_JDay,
                          models = hourly_models,
                          misstolerance = 10)

past_observed <- cbind(
  daily_models_observed,
  hourly_models_observed[[1]][,names(hourly_models)])

save_temperature_scenarios(past_scenarios,
                           "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate", "Wesel_multichill_305_59_historic")
write.csv(past_observed,
          "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate/Wesel_multichill_305_59_observed.csv",
          row.names=FALSE)
SSPs <- c("ssp126", "ssp245","ssp370", "ssp585")
Times <- c(2050, 2085)

future_temps <- load_temperature_scenarios("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate","Wesel_future_")

list_ssp <- 
  strsplit(names(future_temps), '\\.') %>%
  map(2) %>%
  unlist()

list_gcm <-
  strsplit(names(future_temps), '\\.') %>%
  map(3) %>%
  unlist()

list_time <-
  strsplit(names(future_temps), '\\.') %>%
  map(4) %>%
  unlist()
for(SSP in SSPs)
{
  for(Time in Times)
    {
    Temps <- future_temps[list_ssp == SSP & list_time == Time]
    names(Temps) <- list_gcm[list_ssp == SSP & list_time == Time]
    daily_models_future_scenarios <- tempResponse_list_daily(
      Temps,
      Start_JDay = Start_JDay,
      End_JDay = End_JDay,
      models = daily_models)
    daily_models_future_scenarios<-lapply(
      daily_models_future_scenarios,
      function(x) x[which(x$Perc_complete>90),])
    hourly_models_future_scenarios<-
      tempResponse_daily_list(
        Temps,
        latitude = 51.405,
        Start_JDay = Start_JDay,
        End_JDay = End_JDay,
        models=hourly_models,
        misstolerance = 10)

    future_scenarios <- daily_models_future_scenarios
    
    future_scenarios <- lapply(
      names(future_scenarios),
      function(x)
        cbind(future_scenarios[[x]],
              hourly_models_future_scenarios[[x]][,names(hourly_models)]))
    names(future_scenarios)<-names(daily_models_future_scenarios)
    
    chill<-future_scenarios
    
    save_temperature_scenarios(
      chill,
      "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
      paste0("Wesel_multichill_305_59_",Time,"_",SSP))
    }
}
chill_past_scenarios <- load_temperature_scenarios(
  "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
  "Wesel_multichill_305_59_historic")

chill_observed <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate/Wesel_multichill_305_59_observed.csv")


chills <- make_climate_scenario(chill_past_scenarios,
                                caption = "Historical",
                                historic_data = chill_observed,
                                time_series = TRUE)

for(SSP in SSPs)
  for(Time in Times)
    {
    chill <- load_temperature_scenarios(
      "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/future_climate",
      paste0("Wesel_multichill_305_59_",Time,"_",SSP))
    if(SSP == "ssp126") SSPcaption <- "SSP1"
    if(SSP == "ssp245") SSPcaption <- "SSP2"
    if(SSP == "ssp370") SSPcaption <- "SSP3"
    if(SSP == "ssp585") SSPcaption <- "SSP5"    
    if(Time == "2050") Time_caption <- "2050"
    if(Time == "2085") Time_caption <- "2085"
    chills <- make_climate_scenario(chill,
                                    caption = c(SSPcaption,
                                                Time_caption),
                                    add_to = chills)
  }

for(i in 1:length(chills))
   {ch <- chills[[i]]
   if(ch$caption[1] == "Historical")
     {GCMs <- rep("none",length(names(ch$data)))
      SSPs <- rep("none",length(names(ch$data)))
      Years <- as.numeric(ch$labels)
      Scenario <- rep("Historical",
                      length(names(ch$data)))} else
                        {GCMs <- names(ch$data)
                        SSPs <- rep(ch$caption[1],
                                    length(names(ch$data)))
                        Years <- rep(as.numeric(ch$caption[2]),
                                     length(names(ch$data)))
                        Scenario <- rep("Future",
                                        length(names(ch$data)))}
   for(nam in names(ch$data))
     {for(met in metrics)
       {temp_res <-
         data.frame(Metric = met,
                    GCM = GCMs[which(nam == names(ch$data))],
                    SSP = SSPs[which(nam == names(ch$data))],
                    Year = Years[which(nam == names(ch$data))],
                    Result = quantile(ch$data[[nam]][,met],0.1), 
                    Scenario = Scenario[which(nam == names(ch$data))])
       if(i == 1 & nam == names(ch$data)[1] & met == metrics[1])
         results <- temp_res else
           results <- rbind(results,
                            temp_res)
         }
     }
   }

for(met in metrics)
  results[which(results$Metric == met),"SWC"] <-
    results[which(results$Metric == met),"Result"]/
      results[which(results$Metric == met & results$Year == 1990),
              "Result"]-1
  1. Make a heat map illustrating past and future changes in Safe Winter Chill, relative to a past scenario, for the 13 chill models used here.
rng = range(results$SWC)

p_future <- ggplot(results[which(!results$GCM == "none"),],
                   aes(GCM,
                       y = factor(Metric,
                                  levels = metrics),
                       fill = SWC)) +
  geom_tile()

p_future <-
  p_future +
  facet_grid(SSP ~ Year) 

p_future <-
  p_future +
  theme_bw(base_size = 15) +
  theme(axis.text = element_text(size=6))

library(colorRamps)
p_future <-
  p_future +
  scale_fill_gradientn(colours = matlab.like(15),
                       labels = scales::percent,
                       limits = rng)

p_future <-
  p_future  +
  theme(axis.text.x = element_text(angle = 75, 
                                   hjust = 1,
                                   vjust = 1)) +
  labs(fill = "Change in\nSafe Winter Chill\nsince 1980") +
  scale_y_discrete(labels = model_labels) +
  ylab("Chill metric")

p_future

p_past<-
  ggplot(results[which(results$GCM == "none"),],
         aes(Year,
             y = factor(Metric, 
                        levels=metrics),
             fill = SWC)) +
  geom_tile()

p_past<-
  p_past +
  theme_bw(base_size = 15) +
  theme(axis.text = element_text(size = 6))

p_past<-
  p_past +
  scale_fill_gradientn(colours = matlab.like(15),
                       labels = scales::percent,
                       limits = rng)

p_past<-
  p_past +
  scale_x_continuous(position = "top") 

p_past<-
  p_past +
  labs(fill = "Change in\nSafe Winter Chill\nsince 1980") +
  scale_y_discrete(labels = model_labels) +
  ylab("Chill metric")

p_past

chill_comp_plot<-
  (p_past +
     p_future +
     plot_layout(guides = "collect",
                 nrow = 2,
                 heights = c(1,3))) &
  theme(legend.position = "right",
        strip.background = element_blank(),
        strip.text = element_text(face = "bold"))

chill_comp_plot

  1. Produce an animated line plot of your results (summarizing Safe Winter Chill across all the GCMs).
hist_results <- results[which(results$GCM == "none"),]
hist_results$SSP <- "SSP1"
hist_results_2 <- hist_results
hist_results_2$SSP <- "SSP2"
hist_results_3 <- hist_results
hist_results_3$SSP <- "SSP3"
hist_results_4 <- hist_results
hist_results_4$SSP <- "SSP5"
hist_results <- rbind(hist_results,
                      hist_results_2,
                      hist_results_3,
                      hist_results_4)

future_results <- results[which(!results$GCM == "none"),]

GCM_aggregate <- aggregate(
  future_results$SWC,
  by=list(future_results$Metric,
          future_results$SSP,
          future_results$Year),
  FUN=mean)

colnames(GCM_aggregate) <- c("Metric",
                             "SSP",
                             "Year",
                             "SWC")

SSP_Time_series<-rbind(hist_results[,c("Metric",
                                       "SSP",
                                       "Year",
                                       "SWC")],
                       GCM_aggregate)


SSP_Time_series$Year <- as.numeric(SSP_Time_series$Year)

chill_change_plot<-
  ggplot(data = SSP_Time_series,
         aes(x = Year,
             y = SWC,
             col = factor(Metric,
                          levels = metrics))) +
  geom_line(lwd = 1.3) +
  facet_wrap(~SSP,
             nrow = 4) +
  theme_bw(base_size = 10) +
  labs(col = "Change in\nSafe Winter Chill\nsince 1980") +
  scale_color_discrete(labels = model_labels) +
  scale_y_continuous(labels = scales::percent) +
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold")) +
  ylab("Safe Winter Chill")

chill_change_plot

library(gganimate)
library(gifski)
## Warning: Paket 'gifski' wurde unter R Version 4.4.2 erstellt
library(png)
library(transformr)
## Warning: Paket 'transformr' wurde unter R Version 4.4.2 erstellt
ccp<-chill_change_plot +
  transition_reveal(Year)

animate(ccp, fps = 10)
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?

anim_save("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/chill_comparison_animation.gif",
          animation = last_animation())

Chapter 20 Simple phenology analysis

  1. Provide a brief narrative describing what p-hacking is, and why this is a problematic approach to data analysis.

    • p-hacking

      Is a manipulative approach to data analysis aimed at achieving statistically significant results. Common techniques include selective reporting, cherry-picking variables, or running multiple tests until the p-value falls below 0.05.

    • Problems

      It increases false positive results which can create misleading conclusions. It also undermines scientific integrity and makes it impossible or very difficult to reproduce results. Additionally, it encourages to only publish “positive” results. All this leads to wasting resources, it can mislead future research and it lowers the public trust in science.

  2. Provide a sketch of your casual understanding of the relationship between temperature and bloom dates.

    Warmer winter and spring temperatures generally lead to earlier bloom dates. Plants require a certain amount of warmth (growing degree days) to sprout in spring. However, they also need a period of cold (chilling requirement) to properly transition into growth.

    Some species need a specific number of cold days in winter before they can respond to warmth. After meeting the chilling requirement, plants need sufficient warmth to trigger blooming. If one requirement isn’t met the plant may fail to bloom or show reduced productivity.

  3. What do we need to know to build a process-based model from this?

    • Key inputs

      • Temperature data (daily min/max, long-term trends)

      • Species-specific chilling and heat accumulation thresholds

      • Day length (photoperiod) for each day

      • Environmental factors influencing temperature (e.g., topography)

    • Model Structure

      • Accumulate chilling degree days from a defined starting point in autumn

      • After chilling requirements are met, track growing degree days

      • Link daily min/max temperatures with day length to estimate hourly temperature trends

Chapter 21 Simple phenology analysis

  1. Briefly explain why you shouldn’t take the results of a PLS regression analysis between temperature and phenology at face value. What do you need in addition in order to make sense of such outputs?

    PLS is is a useful mining tool, but not a good analytical method on it’s own, especially for small datasets. If there are only limited observations and a high variability, then the results can be misleading when the attempt of an interpretation is made. PLS should therefore always be used with a strong theoretical framework to validate or challenge existing knowledge about e.g. temperature responses.

    Temperature data with a high amount of data (like daily or even hourly data) is a challenge, because the phenology data that it is compared with is very limited (e.g. one bloom date or period per year). Creating a statistical correlation between these two is rather difficult. The effects of temperature vary over time, and simple regression approaches often fail to capture these complexities. That is why PLS results should always be cross validated with other models or experimental data to ensure meaningful interpretations.

  2. Replicate the PLS analysis for the Roter Boskoop dataset that you used in a previous lesson.

Boskoop <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/Roter_Boskoop_bloom_1958_2019.csv") %>%
  select(Pheno_year, First_bloom) %>%
  mutate(Year = as.numeric(substr(First_bloom, 1, 4)),
         Month = as.numeric(substr(First_bloom, 5, 6)),
         Day = as.numeric(substr(First_bloom, 7, 8))) %>%
  make_JDay() %>%
  select(Pheno_year, JDay) %>%
  rename(Year = Pheno_year,
         pheno = JDay)

# head(Boskoop)

KA_temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/TMaxTMin1958-2019_patched.csv") %>%
  make_JDay()

# head(KA_temps)

PLS_results <- PLS_pheno(KA_temps,
                         Boskoop)

# head(PLS_results$PLS_summary)
ggplot_PLS <- function(PLS_results)
{
  library(ggplot2)
  PLS_gg <- PLS_results$PLS_summary %>%
    mutate(Month = trunc(Date / 100),
           Day = Date - Month * 100,
           Date = NULL) 
  
  PLS_gg$Date <- ISOdate(2002, 
                         PLS_gg$Month, 
                         PLS_gg$Day)
  
  PLS_gg$Date[PLS_gg$JDay <= 0] <-
    ISOdate(2001, 
            PLS_gg$Month[PLS_gg$JDay <= 0], 
            PLS_gg$Day[PLS_gg$JDay <= 0])
  
  PLS_gg <- PLS_gg %>%
    mutate(VIP_importance = VIP >= 0.8,
           VIP_Coeff = factor(sign(Coef) * VIP_importance))
  
  VIP_plot<- ggplot(PLS_gg,aes(x=Date,y=VIP)) +
    geom_bar(stat='identity',aes(fill=VIP>0.8))
  
  VIP_plot <- VIP_plot +
    scale_fill_manual(name="VIP", 
                      labels = c("<0.8", ">0.8"), 
                      values = c("FALSE" = "grey", 
                                 "TRUE" = "blue")) +
    theme_bw(base_size=15) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank())
  
  coeff_plot <- ggplot(PLS_gg,
                       aes(x = Date,
                           y = Coef)) +
    geom_bar(stat ='identity',
             aes(fill = VIP_Coeff)) +
    scale_fill_manual(name = "Effect direction", 
                      labels = c("Advancing",
                                 "Unimportant",
                                 "Delaying"), 
                      values = c("-1" = "red", 
                                 "0" = "grey",
                                 "1" = "dark green")) +
    theme_bw(base_size = 15) +
    ylab("PLS coefficient") +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank())
  
  temp_plot <- ggplot(PLS_gg) +
    geom_ribbon(aes(x = Date,
                    ymin = Tmean - Tstdev,
                    ymax = Tmean + Tstdev),
                fill = "grey") +
    geom_ribbon(aes(x = Date,
                    ymin = Tmean - Tstdev * (VIP_Coeff == -1),
                    ymax = Tmean + Tstdev * (VIP_Coeff == -1)),
                fill = "red") +
    geom_ribbon(aes(x = Date,
                    ymin = Tmean - Tstdev * (VIP_Coeff == 1),
                    ymax = Tmean + Tstdev * (VIP_Coeff == 1)),
                fill = "dark green") +
    geom_line(aes(x = Date,
                  y = Tmean)) +
    theme_bw(base_size = 15) +
    ylab(expression(paste(T[mean]," (°C)")))
  
  library(patchwork)
  
  plot<- (VIP_plot +
            coeff_plot +
            temp_plot +
            plot_layout(ncol=1,
                        guides = "collect")
          ) & theme(legend.position = "right",
                    legend.text = element_text(size = 8),
                    legend.title = element_text(size = 10),
                    axis.title.x = element_blank())
  
  plot}

ggplot_PLS(PLS_results)

  1. Write down your thoughts on why we’re not seeing the temperature response pattern we may have expected. What happened to the chill response?

    The expected chilling response is unclear in the PLS analysis, possibly because of some data limitations, overlapping chilling and forcing phases, or the chosen analysis parameters. Chilling accumulates over an extended period of time, and the model may not fully capture this dynamic. Also the smoothing effect of the 11-day running mean or an incorrect start date for the phenological year might also hide some chilling signals.

    Additionally, environmental factors beyond temperature, such as photoperiod or soil moisture, could influence the plants but these are not included as an influence in the model.

Chapter 22 Successes and limitations of PLS regression analysis

  1. Briefly explain in what climatic settings we can expect PLS regression to detect the chilling phase - and in what settings this probably wont work.

    PLS would be able to detect the chilling phase in locations with mild winters where warm temperatures reduce the chill accumulation (like in California). In regions with colder winters (like Beijing, China) where a temperature change can increase and decrease chilling it is way more difficult for PLS to detect chilling phases. This is because if it gets too cold the chill accumulation stops and this creates a dynamic that a simple PLS can not detect.

  2. How could we overcome this problem?

    We could create a model that explicitly calculates the chill accumulation instead of an RLS or we apply alternative statistical methods that can capture the nonlinear relationship between chilling and the temperature.

Chapter 23 PLS regression with agroclimatic metrics

  1. Repeat the PLS_chill_force procedure for the ‘Roter Boskoop’ dataset. Include plots of daily chill and heat accumulation.
temps_hourly <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/TMaxTMin1958-2019_patched.csv") %>%
  stack_hourly_temps(latitude = 50.6)

daychill <- daily_chill(hourtemps = temps_hourly,
                        running_mean = 1,
                        models = list(
                          Chilling_Hours = Chilling_Hours,
                          Utah_Chill_Units = Utah_Model,
                          Chill_Portions = Dynamic_Model,
                          GDH = GDH))

dc <- make_daily_chill_plot2(daychill,
                             metrics = c("Chill_Portions"),
                             cumulative = FALSE,
                             startdate = 300,
                             enddate = 30,
                             focusyears = c(2008), 
                             metriclabels = "Chill Portions")

dc <- make_daily_chill_plot2(daychill,
                             metrics = c("Chill_Portions"),
                             cumulative = TRUE,
                             startdate = 300,
                             enddate = 30,
                             focusyears = c(2008),
                             metriclabels = "Chill Portions")

Boskoop <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/Roter_Boskoop_bloom_1958_2019.csv") %>%
  select(Pheno_year, First_bloom) %>%
  mutate(Year = as.numeric(substr(First_bloom, 1, 4)),
         Month = as.numeric(substr(First_bloom, 5, 6)),
         Day = as.numeric(substr(First_bloom, 7, 8))) %>%
  make_JDay() %>%
  select(Pheno_year, 
         JDay) %>%
  rename(Year = Pheno_year,
         pheno = JDay)

plscf <- PLS_chill_force(daily_chill_obj = daychill,
                         bio_data_frame = Boskoop,
                         split_month = 6,
                         chill_models = "Chill_Portions",
                         heat_models = "GDH")
plot_PLS(plscf,
         PLS_results_path = "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/plscf_outputs")

plscf <- PLS_chill_force(daily_chill_obj = daychill,
                         bio_data_frame = Boskoop,
                         split_month = 6,
                         chill_models = "Chill_Portions",
                         heat_models = "GDH",
                         runn_means = 11)

plot_PLS(plscf,
         PLS_results_path = "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/plscf_outputs_11days")

plot_PLS(plscf,
         PLS_results_path = "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/plscf_outputs_11days_periods",
         add_chill = c(-48,62),
         add_heat = c(3,105.5))

PLS_gg <- plscf$Chill_Portions$GDH$PLS_summary %>%
  mutate(Month = trunc(Date/100),
         Day = Date - Month * 100,
         Date = ISOdate(2002,
                        Month,
                        Day))

PLS_gg[PLS_gg$JDay <= 0,"Date"]<-
  ISOdate(2001,
          PLS_gg$Month[PLS_gg$JDay <= 0],
          PLS_gg$Day[PLS_gg$JDay <= 0])

PLS_gg <- PLS_gg %>%
  mutate(VIP_importance = VIP >= 0.8,
         VIP_Coeff = factor(sign(Coef) * VIP_importance))

chill_start_JDay <- -48
chill_end_JDay <- 62
heat_start_JDay <- 10
heat_end_JDay <- 116.5

chill_start_date <- ISOdate(2001,
                            12,
                            31) + chill_start_JDay * 24 * 3600
chill_end_date <- ISOdate(2001,
                          12,
                          31) + chill_end_JDay * 24 * 3600
heat_start_date <- ISOdate(2001,
                           12,
                           31) + heat_start_JDay * 24 * 3600
heat_end_date <- ISOdate(2001,
                         12,
                         31) + heat_end_JDay * 24 * 3600

temp_plot <- ggplot(PLS_gg,
                    x = Date) +
  annotate("rect",
           xmin = chill_start_date,
           xmax = chill_end_date,
           ymin = -Inf,
           ymax = Inf,
           alpha = .1,
           fill = "blue") +
  annotate("rect",
           xmin = heat_start_date,
           xmax = heat_end_date,
           ymin = -Inf,
           ymax = Inf,
           alpha = .1,
           fill = "red") +
  annotate("rect",
           xmin = ISOdate(2001,
                          12,
                          31) +
             min(plscf$pheno$pheno,
                 na.rm = TRUE) * 24 * 3600,
           xmax = ISOdate(2001,
                          12,
                          31) +
             max(plscf$pheno$pheno,
                 na.rm = TRUE) * 24 * 3600,
           ymin = -Inf,
           ymax = Inf,
           alpha = .1,
           fill = "black") +
  geom_vline(xintercept = ISOdate(2001,
                                  12,
                                  31) +
               median(plscf$pheno$pheno,
                      na.rm = TRUE) * 24 * 3600,
             linetype = "dashed") +
  geom_ribbon(aes(x = Date,
                  ymin = MetricMean - MetricStdev ,
                  ymax = MetricMean + MetricStdev),
              fill="grey") +
  geom_ribbon(aes(x = Date,
                  ymin = MetricMean - MetricStdev * (VIP_Coeff == -1),
                  ymax = MetricMean + MetricStdev * (VIP_Coeff == -1)),
              fill = "red") +
  geom_ribbon(aes(x = Date,
                  ymin = MetricMean - MetricStdev * (VIP_Coeff == 1),
                  ymax = MetricMean + MetricStdev * (VIP_Coeff == 1)),
              fill = "dark green") +
  geom_line(aes(x = Date,
                y = MetricMean ))

temp_plot

plot_PLS_chill_force <- function(plscf,
                                 chill_metric = "Chill_Portions",
                                 heat_metric = "GDH",
                                 chill_label = "CP",
                                 heat_label = "GDH",
                                 chill_phase = c(-48, 62),
                                 heat_phase = c(10, 116.5))
{
  PLS_gg <- plscf[[chill_metric]][[heat_metric]]$PLS_summary %>%
    mutate(Month = trunc(Date/100),
           Day = Date - Month * 100,
           Date = ISOdate(2002,
                          Month,
                          Day))
  
  PLS_gg[PLS_gg$JDay <= 0,"Date"]<-
    ISOdate(2001,
            PLS_gg$Month[PLS_gg$JDay <= 0],
            PLS_gg$Day[PLS_gg$JDay <= 0])
  
  PLS_gg <- PLS_gg %>%
    mutate(VIP_importance = VIP >= 0.8,
           VIP_Coeff = factor(sign(Coef) * VIP_importance))
  
  chill_start_date <- ISOdate(2001,
                              12,
                              31) + chill_phase[1] * 24 * 3600
  chill_end_date <- ISOdate(2001,
                            12,
                            31) + chill_phase[2] * 24 * 3600
  heat_start_date <- ISOdate(2001,
                             12,
                             31) + heat_phase[1] * 24 * 3600
  heat_end_date <- ISOdate(2001,
                           12,
                           31) + heat_phase[2] * 24 * 3600

  temp_plot <- ggplot(PLS_gg) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) +
               min(plscf$pheno$pheno,
                   na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) +
               max(plscf$pheno$pheno,
                   na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) +
                 median(plscf$pheno$pheno,
                        na.rm=TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev ,
                    ymax = MetricMean + MetricStdev ),
                fill = "grey") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev * (VIP_Coeff == -1),
                    ymax = MetricMean + MetricStdev * (VIP_Coeff == -1)),
                fill = "red") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev * (VIP_Coeff == 1),
                    ymax = MetricMean + MetricStdev * (VIP_Coeff == 1)),
                fill = "dark green") +
    geom_line(aes(x = Date,
                  y = MetricMean)) +
    facet_wrap(vars(Type),
               scales = "free_y",
               strip.position = "left",
               labeller = 
                 labeller(
                   Type =
                     as_labeller(c(Chill = paste0("Chill (",
                                                  chill_label,
                                                  ")"),
                                   Heat = paste0("Heat (",
                                                 heat_label,
                                                 ")"))))) +
    ggtitle("Daily chill and heat accumulation rates") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          )
  
  VIP_plot <- ggplot(PLS_gg,
                     aes(x = Date,
                         y = VIP)) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) + min(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) + max(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) + median(plscf$pheno$pheno,
                                                 na.rm = TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_bar(stat = 'identity',
             aes(fill = VIP>0.8)) +
    facet_wrap(vars(Type), 
               scales = "free",
               strip.position = "left",
               labeller = 
                 labeller(
                   Type = as_labeller(c(Chill="VIP for chill",
                                        Heat="VIP for heat")))) +
    scale_y_continuous(
      limits = c(0,
                 max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$VIP))) +
    ggtitle("Variable Importance in the Projection (VIP) scores") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          ) +
    scale_fill_manual(name = "VIP", 
                      labels = c("<0.8", ">0.8"), 
                      values = c("FALSE" = "grey",
                                 "TRUE" = "blue")) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  coeff_plot <- ggplot(PLS_gg,
                       aes(x = Date,
                           y = Coef)) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) + min(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) + max(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) + median(plscf$pheno$pheno,
                                                 na.rm = TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_bar(stat = 'identity',
             aes(fill = VIP_Coeff)) +
    facet_wrap(vars(Type),
               scales = "free",
               strip.position = "left",
               labeller =
                 labeller(
                   Type = as_labeller(c(Chill = "MC for chill",
                                        Heat = "MC for heat")))) +
    scale_y_continuous(
      limits = c(min(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef),
                 max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef))) +
    ggtitle("Model coefficients (MC)") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          ) +
    scale_fill_manual(name = "Effect direction", 
                      labels = c("Advancing",
                                 "Unimportant",
                                 "Delaying"), 
                      values = c("-1" = "red",
                                 "0" = "grey",
                                 "1" = "dark green")) +
    ylab("PLS coefficient") +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  library(patchwork)
  
  plot <- (VIP_plot +
             coeff_plot +
             temp_plot +
             plot_layout(ncol = 1,
                         guides = "collect")
           ) & theme(legend.position = "right",
                     legend.text = element_text(size = 8),
                     legend.title = element_text(size = 10),
                     axis.title.x = element_blank())
plot
}

plot_PLS_chill_force(plscf)

  1. Run PLS_chill_force analyses for all three major chill models. Delineate your best estimates of chilling and forcing phases for all of them.

  2. Plot results for all three analyses, including shaded plot areas for the chilling and forcing periods you estimated.

daychill <- daily_chill(hourtemps = temps_hourly,
                        running_mean = 11,
                        models = list(Chilling_Hours = Chilling_Hours,
                                      Utah_Chill_Units = Utah_Model,
                                      Chill_Portions = Dynamic_Model,
                                      GDH = GDH))

plscf <- PLS_chill_force(daily_chill_obj = daychill,
                         bio_data_frame = Boskoop,
                         split_month = 6,
                         chill_models = c("Chilling_Hours",
                                          "Utah_Chill_Units",
                                          "Chill_Portions"),
                         heat_models = c("GDH"))

plot_PLS_chill_force(plscf,
                     chill_metric = "Chilling_Hours",
                     heat_metric = "GDH",
                     chill_label = "CH",
                     heat_label = "GDH",
                     chill_phase = c(-30, 13),
                     heat_phase = c(-10, 116.5))

plot_PLS_chill_force(plscf,
                     chill_metric = "Utah_Chill_Units",
                     heat_metric = "GDH",
                     chill_label = "CU",
                     heat_label = "GDH",
                     chill_phase = c(-70, 20),
                     heat_phase = c(10, 116.5))

plot_PLS_chill_force(plscf,
                     chill_metric = "Chill_Portions",
                     heat_metric = "GDH",
                     chill_label = "CP",
                     heat_label = "GDH",
                     chill_phase = c(-75, 20),
                     heat_phase = c(0, 116.5))

Chapter 24 Examples of PLS regression with agroclimatic metrics

  1. Look across all the PLS results presented above. Can you detect a pattern in where chilling and forcing periods delineated clearly, and where this attempt failed?

    • Chilling period

      • In warm climates like in Croatia or Tunisia is the chilling period is clearly visible.

        • The chilling period is only as long as it needs to be or sometimes even too short what makes it a lot easier to detect the chilling period.
      • In cold climates like in Beijing or even Klein-Altendorf is the chilling period not easy to detect.

        • The chilling period is longer than required, what makes it difficult to detect a clear chilling phase
    • Forcing period

      • It is easy to detect the forcing period in basically all climates.

      • This is probably do to the fact that the heat accumulation has a more direct influence on blooming than chilling.

  2. Think about possible reasons for the success or failure of PLS based on agroclimatic metrics. Write down your thoughts.

    • Reasons for success

      • Detecting the two phases is easier in climates that have a clearer separation between the chilling and the forcing phase.

      • If the chilling requirements are just met then the response is clearer to detect.

      • The forcing period is nearly always good to detect because it has a more direct influence to blooming than chilling.

      • The choice of Chill model also has an influence on the results, a more simple model like Chilling Hours get less satisfying results than the Dynamic Model.

    • Reasons for failure

      • A rather cold climate with a longer chilling period and an overlapping chilling and forcing period is difficult to analyse with a PLS analysis.

      • PLS analysis assumes linear relationships, which may not capture the complex interactions between chilling, heat, and bloom timing. If only forcing period and blooming are viewed there is a more linear relationship that PLS can detect.

Chapter 25 Why PLS doesn’t always work

  1. Produce chill and heat model sensitivity plots for the location you focused on in previous exercises.
mon <- 1
ndays <- 31
tmin <- 1
tmax <- 8
latitude <- 51.4


weather <- make_all_day_table(
            data.frame(Year = c(2001, 2001),
                       Month = c(mon, mon),
                       Day = c(1, ndays),
                       Tmin = c(0, 0),
                       Tmax = c(0, 0))) %>%
  mutate(Tmin = tmin,
         Tmax = tmax)

hourly_temps <- stack_hourly_temps(weather,
                                   latitude = latitude)

CPs <- Dynamic_Model(
     hourly_temps$hourtemps$Temp)

daily_CPs <- CPs[length(CPs)] / nrow(weather)

# daily_CPs


month_range <- c(10, 11, 12, 1, 2, 3)

Tmins <- c(-20:20)
Tmaxs <- c(-15:30)

mins <- NA
maxs <- NA
CP <- NA
month <- NA
temp_model <- Dynamic_Model

for(mon in month_range)
    {days_month <- as.numeric(
      difftime( ISOdate(2002,
                        mon + 1,
                        1),
                ISOdate(2002,
                        mon,
                        1)))
    
    if(mon == 12) days_month <- 31
    
    weather <- make_all_day_table(
                data.frame(Year = c(2001, 2001),
                           Month = c(mon, mon),
                           Day = c(1, days_month),
                           Tmin = c(0, 0),
                           Tmax = c(0, 0)))
    
    for(tmin in Tmins)
      for(tmax in Tmaxs)
        if(tmax >= tmin)
          {
          hourtemps <- weather %>%
            mutate(Tmin = tmin,
                   Tmax = tmax) %>%
            stack_hourly_temps(latitude = latitude) %>% 
            pluck("hourtemps", "Temp")

          CP <- c(CP,
                  tail(do.call(temp_model,
                               list(hourtemps)), 1) /
                    days_month)
          mins <- c(mins, tmin)
          maxs <- c(maxs, tmax)
          month <- c(month, mon)
        }
}

results <- data.frame(Month = month,
                      Tmin = mins,
                      Tmax = maxs,
                      CP)

results <- results[!is.na(results$Month), ]


write.csv(results,
          "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/model_sensitivity_development.csv",
          row.names = FALSE)

results$Month_names <- factor(results$Month,
                              levels = month_range,
                              labels = month.name[month_range])  

DM_sensitivity <- ggplot(results,
                         aes(x = Tmin,
                             y = Tmax,
                             fill = CP)) +
  geom_tile() +
  scale_fill_gradientn(colours = alpha(matlab.like(15),
                                       alpha = .5),
                       name = "Chill/day (CP)") +
  ylim(min(results$Tmax),
       max(results$Tmax)) +
  ylim(min(results$Tmin),
       max(results$Tmin))

DM_sensitivity

DM_sensitivity <- DM_sensitivity +
  facet_wrap(vars(Month_names)) +
  ylim(min(results$Tmax),
       max(results$Tmax)) +
  ylim(min(results$Tmin),
       max(results$Tmin))

DM_sensitivity

latitude <- 51.4

month_range <- c(10, 11, 12, 1, 2, 3)

temperatures <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/patched_weather_3_Wesel.csv") %>% 
  filter(Month %in% month_range) %>%
  mutate(Month_names =
           factor(Month,
                  levels = c(10, 11, 12, 1, 2, 3),
                  labels = c("October", "November", "December",
                             "January", "February", "March")))

temperatures[which(temperatures$Tmax < temperatures$Tmin),
             c("Tmax", "Tmin")] <- NA

DM_sensitivity +
  geom_point(data = temperatures,
             aes(x = Tmin,
                 y = Tmax,
                 fill = NULL,
                 color = "Temperature"),
             size = 0.2) +
  facet_wrap(vars(Month_names)) +
  scale_color_manual(values = "black",
                     labels = "Daily temperature \nextremes (°C)",
                     name = "Observed at site" ) +
  guides(fill = guide_colorbar(order = 1),
         color = guide_legend(order = 2)) +
  ylab("Tmax (°C)") +
  xlab("Tmin (°C)") + 
  theme_bw(base_size = 15) 

Chapter 26 Evaluating PLS outputs

  1. Reproduce the analysis for the ‘Roter Boskoop’ dataset.
roter_B <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/Roter_Boskoop_bloom_1958_2019.csv") %>% 
  select(Pheno_year, First_bloom) %>%
  mutate(Year = as.numeric(substr(First_bloom, 1, 4)),
         Month = as.numeric(substr(First_bloom, 5, 6)),
         Day = as.numeric(substr(First_bloom, 7, 8))) %>%
  make_JDay() %>%
  select(Pheno_year, 
         JDay) %>%
  rename(Year = Pheno_year,
         pheno = JDay)

temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/TMaxTMin1958-2019_patched.csv")

temps_hourly <- temps %>%
  stack_hourly_temps(latitude = 50.6)

daychill <- daily_chill(hourtemps = temps_hourly,
                        running_mean = 1,
                        models = list(Chilling_Hours = Chilling_Hours,
                                      Utah_Chill_Units = Utah_Model,
                                      Chill_Portions = Dynamic_Model,
                                      GDH = GDH)
                        )


plscf <- PLS_chill_force(daily_chill_obj = daychill,
                         bio_data_frame = roter_B,
                         split_month = 6,
                         chill_models = "Chill_Portions",
                         heat_models = "GDH",
                         runn_means = 11)
plot_PLS_chill_force <- function(plscf,
                                 chill_metric = "Chill_Portions",
                                 heat_metric = "GDH",
                                 chill_label = "CP",
                                 heat_label = "GDH",
                                 chill_phase = c(-48, 62),
                                 heat_phase = c(10, 116.5))
{
  PLS_gg <- plscf[[chill_metric]][[heat_metric]]$PLS_summary %>%
    mutate(Month = trunc(Date/100),
           Day = Date - Month * 100,
           Date = ISOdate(2002,
                          Month,
                          Day))
  
  PLS_gg[PLS_gg$JDay <= 0,"Date"]<-
    ISOdate(2001,
            PLS_gg$Month[PLS_gg$JDay <= 0],
            PLS_gg$Day[PLS_gg$JDay <= 0])
  
  PLS_gg <- PLS_gg %>%
    mutate(VIP_importance = VIP >= 0.8,
           VIP_Coeff = factor(sign(Coef) * VIP_importance))
  
  chill_start_date <- ISOdate(2001,
                              12,
                              31) + chill_phase[1] * 24 * 3600
  chill_end_date <- ISOdate(2001,
                            12,
                            31) + chill_phase[2] * 24 * 3600
  heat_start_date <- ISOdate(2001,
                             12,
                             31) + heat_phase[1] * 24 * 3600
  heat_end_date <- ISOdate(2001,
                           12,
                           31) + heat_phase[2] * 24 * 3600

  temp_plot <- ggplot(PLS_gg) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) +
               min(plscf$pheno$pheno,
                   na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) +
               max(plscf$pheno$pheno,
                   na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) +
                 median(plscf$pheno$pheno,
                        na.rm=TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev ,
                    ymax = MetricMean + MetricStdev ),
                fill = "grey") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev * (VIP_Coeff == -1),
                    ymax = MetricMean + MetricStdev * (VIP_Coeff == -1)),
                fill = "red") +
    geom_ribbon(aes(x = Date,
                    ymin = MetricMean - MetricStdev * (VIP_Coeff == 1),
                    ymax = MetricMean + MetricStdev * (VIP_Coeff == 1)),
                fill = "dark green") +
    geom_line(aes(x = Date,
                  y = MetricMean)) +
    facet_wrap(vars(Type),
               scales = "free_y",
               strip.position = "left",
               labeller = 
                 labeller(
                   Type =
                     as_labeller(c(Chill = paste0("Chill (",
                                                  chill_label,
                                                  ")"),
                                   Heat = paste0("Heat (",
                                                 heat_label,
                                                 ")"))))) +
    ggtitle("Daily chill and heat accumulation rates") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          )
  
  VIP_plot <- ggplot(PLS_gg,
                     aes(x = Date,
                         y = VIP)) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) + min(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) + max(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) + median(plscf$pheno$pheno,
                                                 na.rm = TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_bar(stat = 'identity',
             aes(fill = VIP>0.8)) +
    facet_wrap(vars(Type), 
               scales = "free",
               strip.position = "left",
               labeller = 
                 labeller(
                   Type = as_labeller(c(Chill="VIP for chill",
                                        Heat="VIP for heat")))) +
    scale_y_continuous(
      limits = c(0,
                 max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$VIP))) +
    ggtitle("Variable Importance in the Projection (VIP) scores") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          ) +
    scale_fill_manual(name = "VIP", 
                      labels = c("<0.8", ">0.8"), 
                      values = c("FALSE" = "grey",
                                 "TRUE" = "blue")) +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  coeff_plot <- ggplot(PLS_gg,
                       aes(x = Date,
                           y = Coef)) +
    annotate("rect",
             xmin = chill_start_date,
             xmax = chill_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "blue") +
    annotate("rect",
             xmin = heat_start_date,
             xmax = heat_end_date,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "red") +
    annotate("rect",
             xmin = ISOdate(2001,
                            12,
                            31) + min(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             xmax = ISOdate(2001,
                            12,
                            31) + max(plscf$pheno$pheno,
                                      na.rm = TRUE) * 24 * 3600,
             ymin = -Inf,
             ymax = Inf,
             alpha = .1,
             fill = "black") +
    geom_vline(xintercept = ISOdate(2001,
                                    12,
                                    31) + median(plscf$pheno$pheno,
                                                 na.rm = TRUE) * 24 * 3600,
               linetype = "dashed") +
    geom_bar(stat = 'identity',
             aes(fill = VIP_Coeff)) +
    facet_wrap(vars(Type),
               scales = "free",
               strip.position = "left",
               labeller =
                 labeller(
                   Type = as_labeller(c(Chill = "MC for chill",
                                        Heat = "MC for heat")))) +
    scale_y_continuous(
      limits = c(min(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef),
                 max(plscf[[chill_metric]][[heat_metric]]$PLS_summary$Coef))) +
    ggtitle("Model coefficients (MC)") +
    theme_bw(base_size = 15) + 
    theme(strip.background = element_blank(),
          strip.placement = "outside",
          strip.text.y = element_text(size = 12),
          plot.title = element_text(hjust = 0.5),
          axis.title.y = element_blank()
          ) +
    scale_fill_manual(name = "Effect direction", 
                      labels = c("Advancing",
                                 "Unimportant",
                                 "Delaying"), 
                      values = c("-1" = "red",
                                 "0" = "grey",
                                 "1" = "dark green")) +
    ylab("PLS coefficient") +
    theme(axis.text.x = element_blank(),
          axis.ticks.x = element_blank(),
          axis.title.x = element_blank(),
          axis.title.y = element_blank())
  
  library(patchwork)
  
  plot <- (VIP_plot +
             coeff_plot +
             temp_plot +
             plot_layout(ncol = 1,
                         guides = "collect")
           ) & theme(legend.position = "right",
                     legend.text = element_text(size = 8),
                     legend.title = element_text(size = 10),
                     axis.title.x = element_blank())
plot
}
plot_PLS_chill_force(plscf,
                     chill_metric = "Chill_Portions",
                     heat_metric = "GDH",
                     chill_label = "CP",
                     heat_label = "GDH",
                     chill_phase = c(-48, 62),
                     heat_phase = c(3, 105.5))

chill_phase <- c(317, 62)
heat_phase <- c(3, 105.5)

chill <- tempResponse(hourtemps = temps_hourly,
                      Start_JDay = chill_phase[1],
                      End_JDay = chill_phase[2],
                      models = list(Chill_Portions = Dynamic_Model),
                      misstolerance = 10)

heat <- tempResponse(hourtemps = temps_hourly,
                     Start_JDay = heat_phase[1],
                     End_JDay = heat_phase[2],
                     models = list(GDH = GDH))

ggplot(data = chill,
       aes(x = Chill_Portions)) +
  geom_histogram() +
  ggtitle("Chill accumulation during endodormancy (Chill Portions)") +
  xlab("Chill accumulation (Chill Portions)") +
  ylab("Frequency between 1958 and 2019") +
  theme_bw(base_size = 12)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggplot(data = heat,
       aes(x = GDH)) +
  geom_histogram() +
  ggtitle("Heat accumulation during ecodormancy (GDH)") +
  xlab("Heat accumulation (Growing Degree Hours)") +
  ylab("Frequency between 1958 and 2019") +
  theme_bw(base_size = 12)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

chill_requirement <- mean(chill$Chill_Portions)
chill_req_error <- sd(chill$Chill_Portions)

heat_requirement <- mean(heat$GDH)
heat_req_error <- sd(heat$GDH)
mean_temp_period <- function(
    temps,
    start_JDay,
    end_JDay,
    end_season = end_JDay)
{ temps_JDay <- make_JDay(temps) %>%
  mutate(Season =Year)
  
  if(start_JDay > end_season)
    temps_JDay$Season[which(temps_JDay$JDay >= start_JDay)]<-
        temps_JDay$Year[which(temps_JDay$JDay >= start_JDay)]+1
  
  if(start_JDay > end_season)
    sub_temps <- subset(temps_JDay,
                        JDay <= end_JDay | JDay >= start_JDay)
  
  if(start_JDay <= end_JDay)
    sub_temps <- subset(temps_JDay,
                        JDay <= end_JDay & JDay >= start_JDay)
  
  mean_temps <- aggregate(sub_temps[, c("Tmin", "Tmax")],
                          by = list(sub_temps$Season),
                          FUN = function(x) mean(x,
                                                 na.rm=TRUE))
  
  mean_temps[, "n_days"] <- aggregate(sub_temps[, "Tmin"],
                                      by = list(sub_temps$Season),
                                      FUN = length)[,2]
  
  mean_temps[, "Tmean"] <- (mean_temps$Tmin + mean_temps$Tmax) / 2
  mean_temps <- mean_temps[, c(1, 4, 2, 3, 5)]
  colnames(mean_temps)[1] <- "End_year"
  
  return(mean_temps)
}

mean_temp_chill <- mean_temp_period(temps = temps,
                                    start_JDay = chill_phase[1],
                                    end_JDay = chill_phase[2],
                                    end_season = 60)

mean_temp_heat <- mean_temp_period(temps = temps,
                                   start_JDay = heat_phase[1],
                                   end_JDay = heat_phase[2],
                                   end_season = 60)
mean_temp_chill <- 
  mean_temp_chill[which(mean_temp_chill$n_days >=
                          max(mean_temp_chill$n_days)-1),]
mean_temp_heat <-
  mean_temp_heat[which(mean_temp_heat$n_days >=
                         max(mean_temp_heat$n_days)-1),]

mean_chill <- mean_temp_chill[, c("End_year",
                                  "Tmean")]
colnames(mean_chill)[2] <- "Tmean_chill"

mean_heat <- mean_temp_heat[,c("End_year",
                               "Tmean")]
colnames(mean_heat)[2] <- "Tmean_heat"

phase_Tmeans <- merge(mean_chill,
                      mean_heat, 
                      by = "End_year")


pheno <- roter_B
colnames(pheno)[1] <- "End_year"

Tmeans_pheno <- merge(phase_Tmeans,
                      pheno,
                      by = "End_year")
head(Tmeans_pheno)
##   End_year Tmean_chill Tmean_heat pheno
## 1     1959    2.924324  4.5815534   104
## 2     1960    3.147206  4.5100953   115
## 3     1961    4.047768  6.0223301    99
## 4     1962    2.541892  2.6927184   134
## 5     1963   -3.529730 -0.9257282   127
## 6     1964    1.443243  2.4723301   126
library(fields)
## Lade nötiges Paket: spam
## Spam version 2.11-0 (2024-10-03) is loaded.
## Type 'help( Spam)' or 'demo( spam)' for a short introduction 
## and overview of this package.
## Help for individual functions is also obtained by adding the
## suffix '.spam' to the function name, e.g. 'help( chol.spam)'.
## 
## Attache Paket: 'spam'
## Die folgenden Objekte sind maskiert von 'package:base':
## 
##     backsolve, forwardsolve
## Lade nötiges Paket: viridisLite
## 
## Try help(fields) to get started.
k <- Krig(x = as.matrix(
                Tmeans_pheno[,
                             c("Tmean_chill",
                               "Tmean_heat")]),
          Y = Tmeans_pheno$pheno)

pred <- predictSurface(k)
colnames(pred$z) <- pred$y
rownames(pred$z) <- pred$x

library(reshape2)
## 
## Attache Paket: 'reshape2'
## Das folgende Objekt ist maskiert 'package:tidyr':
## 
##     smiths
melted <- melt(pred$z)
  
library(metR)
## Warning: Paket 'metR' wurde unter R Version 4.4.2 erstellt
## 
## Attache Paket: 'metR'
## Das folgende Objekt ist maskiert 'package:purrr':
## 
##     cross
library(colorRamps)
  
colnames(melted) <- c("Tmean_chill",
                      "Tmean_heat",
                      "value")


ggplot(melted,
       aes(x = Tmean_chill,
           y = Tmean_heat,
           z = value)) +
  geom_contour_fill(bins = 100) +
  scale_fill_gradientn(colours = alpha(matlab.like(15)),
                       name = "Bloom date \n(day of the year)") +
  geom_contour(col = "black")  +
  geom_point(data = Tmeans_pheno,
             aes(x = Tmean_chill,
                 y = Tmean_heat,
                 z = NULL),
             size = 0.7) +
  geom_text_contour(stroke = 0.2) +
  ylab(expression(paste("Forcing phase ", 
                        T[mean],
                        " (",
                        degree,
                        "C)"))) +
  xlab(expression(paste("Chilling phase ",
                        T[mean],
                        " (",
                        degree,
                        "C)")))  +
  theme_bw(base_size = 15)
## Warning: Removed 4368 rows containing non-finite outside the scale range
## (`stat_contour()`).

temps <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/TMaxTMin1958-2019_patched.csv")

roter_B <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/data ch20-25/Roter_Boskoop_bloom_1958_2019.csv")  %>% 
  select(Pheno_year, First_bloom) %>%
  mutate(Year = as.numeric(substr(First_bloom, 1, 4)),
         Month = as.numeric(substr(First_bloom, 5, 6)),
         Day = as.numeric(substr(First_bloom, 7, 8))) %>%
  make_JDay() %>%
  select(Pheno_year, 
         JDay) %>%
  rename(Year = Pheno_year,
         pheno = JDay)

colnames(roter_B) <- c("Year", "pheno")

Cali_temps_hourly <- stack_hourly_temps(temps,
                                        latitude = 50.6)

Cali_daychill <- daily_chill(hourtemps = Cali_temps_hourly,
                             running_mean = 1,
                             models = list(Chilling_Hours = Chilling_Hours,
                                           Utah_Chill_Units = Utah_Model,
                                           Chill_Portions = Dynamic_Model,
                                           GDH = GDH)
    )


plscf <- PLS_chill_force(daily_chill_obj = Cali_daychill,
                         bio_data_frame = roter_B,
                         split_month = 6,
                         chill_models = "Chill_Portions",
                         heat_models = "GDH",
                         runn_means = 11)

plot_PLS_chill_force(plscf,
                     chill_metric = "Chill_Portions",
                     heat_metric = "GDH",
                     chill_label = "CP",
                     heat_label = "GDH",
                     chill_phase = c(-56, 5),
                     heat_phase = c(19, 77))

  1. We’ve looked at data from a number of locations so far. How would you expect this surface plot to look like in Beijing? And how should it look in Tunisia?
    • Beijing
      • The chilling period (x-axis) will have much colder temperatures than Klein-Altendorf.
      • The forcing period (y-axis) will be more narrow than Klein-Altendorf.
      • The separation lines for the bloom dates will be horizontal.
      • The forcing period will mainly influence the bloom date.
    • Tunisia
      • The chilling period (x-axis) will have warmer temperatures than Klein-Altendorf.
      • The forcing period (y-axis) will be a lot wider than Klein-Altendorf.
      • The separation lines for the bloom dates will be nearly vertical.
      • The chilling period will mainly influence the bloom date.

Chapter 27The relative importance of chill and heat

  1. Describe the temperature response hypothesis outlined in this chapter.

    The temperature response hypothesis explains how temperature influences tree phenology through chilling and forcing effects. Chilling temperatures are required to release bud dormancy and forcing temperatures to drive bud development and growth. The hypothesis examines how the balance between chilling and forcing temperatures determines the timing of phenological events. Tree species and even single varieties have specific needs in chilling and forcing to produce sufficient fruits.

    In the graphic below is a representation how the relationship of chilling and forcing influences the bloom date of trees.

Graphical representation of the temperature response hypothesis
Graphical representation of the temperature response hypothesis

Chapter 28 Experimentally enhanced PLS

There are no tasks for this chapter.

Chapter 29

  1. Explain the difference between output validation and process validation.
  • Output validation

    • Compares the model’s outputs to real-world observations to assess accuracy

    • The focus is on the end results that the model produces

  • Process validation

    • Looks into the internal mechanisms and assumptions of the model

    • The focus is on the model’s processes and if they correctly represent real-world systems.

    • So rather than validating only the model’s results it checks if the path to get the results aligns with reality

  1. Explain what a validity domain is and why it is important to consider this whenever we want to use our model to forecast something.
  • Validity domain

    Defines the boundaries within which the model works in a useful way and delivers reliable results.

  • Importance

    Ensures that the model is applied only to situations where it has been validated and prevents misuse.

  1. What is validation for purpose?

    Assesses if a model is suitable for a specific intended application. It checks if the model meets the requirements and objectives of the intended use. Also the accuracy of the model and its relevance to the specific context are evaluated.

  2. How can we ensure that our model is suitable for the predictions we want to make?

    To ensure this there are some steps that need to be followed:

    1. Define the models purpose and clearly outline the specific objectives and applications of the desired model.

    2. Validate your chosen model. For this use both output and process validation to assess accuracy and internal consistency.

    3. Determine the model’s valid range. Identify the conditions under which the model operates reliably.

    4. Do continuous evaluations of the model to ensure its accuracy. Compare the model’s predictions to real-world data and update the model if errors occur.

Chapter 30 The PhenoFlex Model

  1. Parameterize the PhenoFlex model for `Roter Boskoop’ apples.
Boskoop <-
  read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Roter_Boskoop_bloom_1958_2019.csv") %>%
  select(Pheno_year, First_bloom) %>%
  mutate(Year = as.numeric(substr(First_bloom, 1, 4)),
         Month = as.numeric(substr(First_bloom, 5, 6)),
         Day = as.numeric(substr(First_bloom, 7, 8))) %>%
  make_JDay() %>%
  select(Pheno_year, JDay) %>%
  rename(Year = Pheno_year,
         pheno = JDay)

hourtemps <- 
  read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/TMaxTMin1958-2019_patched.csv") %>%
  stack_hourly_temps(latitude = 50.6)

par <-   c(40, 190, 0.5, 25, 3372.8,  9900.3, 6319.5, 5.939917e13,  4, 36,  4,  1.60)
upper <- c(41, 200, 1.0, 30, 4000.0, 10000.0, 7000.0, 6.e13, 10, 40, 10, 50.00)
lower <- c(38, 180, 0.1, 0 , 3000.0,  9000.0, 6000.0, 5.e13,  0,  0,  0,  0.05)

SeasonList <- genSeasonList(hourtemps$hourtemps,
                            mrange = c(8, 6),
                            years = c(1959:2019))
Fit_res <- phenologyFitter(par.guess = par, 
                           modelfn = PhenoFlex_GDHwrapper,
                           bloomJDays = Boskoop$pheno[which(Boskoop$Year > 1957)],
                           SeasonList = SeasonList,
                           lower = lower,
                           upper = upper,
                           control = list(smooth = FALSE,
                                          verbose = FALSE, 
                                          maxit = 1000,
                                          nb.stop.improvement = 5))

Boskoop_par <- Fit_res$par

write.csv(Boskoop_par, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/PhenoFlex_parameters_roter_Boskoop.csv")
  1. Produce plots of predicted vs. observed bloom dates and distribution of prediction errors.
Boskoop_par <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/PhenoFlex_parameters_roter_Boskoop.csv")[,2]

SeasonList <- genSeasonList(hourtemps$hourtemps, 
                            mrange = c(8, 6),
                            years = c(1959:2019))

Boskoop_PhenoFlex_predictions <- Boskoop[which(Boskoop$Year > 1958),]

for(y in 1:length(Boskoop_PhenoFlex_predictions$Year))
   Boskoop_PhenoFlex_predictions$predicted[y] <-
    PhenoFlex_GDHwrapper(SeasonList[[y]],
                         Boskoop_par)

Boskoop_PhenoFlex_predictions$Error <- 
  Boskoop_PhenoFlex_predictions$predicted - 
  Boskoop_PhenoFlex_predictions$pheno

ggplot(Boskoop_PhenoFlex_predictions,
       aes(x = pheno,
           y = predicted)) +
  geom_point() +
  geom_abline(intercept = 0,
              slope = 1) +
  theme_bw(base_size = 15) +
  xlab("Observed bloom date (Day of the year)") +
  ylab("Predicted bloom date (Day of the year)") +
  ggtitle("Predicted vs. observed bloom dates")

ggplot(Boskoop_PhenoFlex_predictions,
       aes(Error)) +
  geom_histogram() +
  ggtitle("Distribution of prediction errors")
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

  1. Compute the model performance metrics RMSEP, mean error and mean absolute error.
RMSEP(Boskoop_PhenoFlex_predictions$predicted,
      Boskoop_PhenoFlex_predictions$pheno)
## [1] 4.52457
mean(Boskoop_PhenoFlex_predictions$Error)
## [1] -1.23125
mean(abs(Boskoop_PhenoFlex_predictions$Error))
## [1] 3.347917

Chapter 31 The PhenoFlex Model - a second look

  1. Make chill and heat response plots for the ‘Roter Boskoop’ PhenoFlex model for the location you did the earlier analyses for.

I really tried to incorporate my own temperature data, but I wasn’t able to fin out how to fix the code to the shorter time span (1990-2019 instead of 1958-2019). Even if I got the bloomJDays andSeasonList to the same length I still got error messages that there are some statements not fully met. So Iswitched back to the CKA temperatures to be able to create the plots correctly.

Boskoop <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Roter_Boskoop_bloom_1958_2019.csv") %>%
  select(Pheno_year, First_bloom) %>%
  mutate(Year = as.numeric(substr(First_bloom, 1, 4)),
         Month = as.numeric(substr(First_bloom, 5, 6)),
         Day = as.numeric(substr(First_bloom, 7, 8))) %>%
  make_JDay() %>%
  select(Pheno_year, JDay) %>%
  rename(Year = Pheno_year,
         pheno = JDay)

temp_data <- read.csv("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/TMaxTMin1958-2019_patched.csv") %>% stack_hourly_temps(latitude = 50.6)

par <-   c(40, 190, 0.5, 25, 3372.8,  9900.3, 6319.5, 5.939917e13,  4, 36,  4,  1.60)
upper <- c(41, 200, 1.0, 30, 4000.0, 10000.0, 7000.0, 6.e13, 10, 40, 10, 50.00)
lower <- c(38, 180, 0.1, 0 , 3000.0,  9000.0, 6000.0, 5.e13,  0,  0,  0,  0.05)

SeasonList <- genSeasonList(temp_data$temp_data,
                            mrange = c(8, 6),
                            years = c(1959:2019))

Fit_res <- phenologyFitter(par.guess = par, 
                           modelfn = PhenoFlex_GDHwrapper,
                           bloomJDays = Boskoop$pheno[which(Boskoop$Year > 1957)],
                           SeasonList = SeasonList,
                           lower = lower,
                           upper = upper,
                           control = list(smooth = FALSE,
                                          verbose = FALSE, 
                                          maxit = 1000,
                                          nb.stop.improvement = 5))

Boskoop_par <- Fit_res$par

write.csv(Boskoop_par, "D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/PhenoFlex_parameters_roter_Boskoop.csv")
apply_const_temp <-
  function(temp, A0, A1, E0, E1, Tf, slope, portions = 1200, deg_celsius = TRUE)
    {
  temp_vector <- rep(temp,
                     times = portions)
  res <- chillR::DynModel_driver(temp = temp_vector,
                                 A0 = A0,
                                 A1 = A1,
                                 E0 = E0,
                                 E1 = E1,
                                 Tf = Tf,
                                 slope = slope,
                                 deg_celsius = deg_celsius)
  return(res$y[length(res$y)])
}

gen_bell <- function(par,
                     temp_values = seq(-5, 20, 0.1)) {
  E0 <- par[5]
  E1 <- par[6]
  A0 <- par[7]
  A1 <- par[8]
  Tf <- par[9]
  slope <- par[12]

  y <- c()
  for(i in seq_along(temp_values)) {
    y[i] <- apply_const_temp(temp = temp_values[i],
                             A0 = A0,
                             A1 = A1,
                             E0 = E0,
                             E1 = E1,
                             Tf = Tf,
                             slope = slope)
  }
  return(invisible(y))
}

GDH_response <- function(T, par)
  {Tb <- par[11]
   Tu <- par[4]
   Tc <- par[10]
   GDH_weight <- rep(0, length(T))
   GDH_weight[which(T >= Tb & T <= Tu)] <-
     1/2 * (1 + cos(pi + pi * (T[which(T >= Tb & T <= Tu)] - Tb)/(Tu - Tb)))
   GDH_weight[which(T > Tu & T <= Tc)] <-
     (1 + cos(pi/2 + pi/2 * (T[which(T >  Tu & T <= Tc)] -Tu)/(Tc - Tu)))
  return(GDH_weight)
}
Boskoop <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/PhenoFlex_parameters_roter_Boskoop.csv")[,2]

temp_values = seq(-5, 30, 0.1)

temp_response <- data.frame(Temperature = temp_values,
                            Chill_response = gen_bell(Boskoop,
                                                      temp_values),
                            Heat_response = GDH_response(temp_values,
                                                         Boskoop))

pivoted_response <- pivot_longer(temp_response, 
                                 c(Chill_response,
                                   Heat_response))

ggplot(pivoted_response,
       aes(x = Temperature,
           y = value)) +
  geom_line(linewidth = 2,
            aes(col = name)) +
  ylab("Temperature response (arbitrary units)") +
  xlab("Temperature (°C)") +
  facet_wrap(vars(name),
             scales = "free",
             labeller = 
               labeller(name = c(Chill_response = c("Chill response"),
                                     Heat_response = c("Heat response")))) +
  scale_color_manual(values = c("Chill_response" = "blue",
                                "Heat_response" = "red")) +
  theme_bw(base_size = 15) +
  theme(legend.position = "none")

Chapter 32 Can we improve the performance of PhenoFlex?

  1. What was the objective of this work?

    The study aimed to enhance the performance of the PhenoFlex model by refining its calibration process.

  2. What was the main conclusion?

    That the PhenoFlex model’s performance improved when using a calibration approach that incorporated the standard deviation and the mean absolute error of phenological dates.

  3. What experiments could we conduct to test the hypothesis that emerged at the end of the conclusion?

    • Cross-validation

      • Apply the model to diverse datasets (different geographic regions and spiecies) to test the robustness and generalizability of the PhenoFlex Model.

      • Use the new data to test the calibrations that were done with the already present data.

    • Comparative analysis

      • Compare the predictive accuracy of the PhenoFlex model calibrated with the combined approach against models calibrated with only one of the error metrics to evaluate the performance differences.
    • Sensitivity analysis

      • Investigate how variations in weighting between standard deviation and mean absolute error during calibration affect the model’s predictions to determine optimal calibration settings.
    • Longitudinal studies

      • Monitor the model’s predictive performance over multiple years to ensure that the improved calibration method maintains its accuracy over time.

Chapter 33 Frost risk analysis

  1. Illustrate the development of the bloom period over the duration of the weather record. Use multiple ways to show this - feel free to be creative.
CKA_Boskoop <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/Roter_Boskoop_bloom_1958_2019.csv")

CKA_weather <- read_tab("D:/! Ricardas Daten/Studium/Master/03. Semester (WiSe 24-25)/Tree phenology analysis in R/Logbook chapters/TMaxTMin1958-2019_patched.csv")

CKA_Boskoop <- 
  CKA_Boskoop %>%
  pivot_longer(cols = "First_bloom":"Last_bloom",
               names_to = "variable",
               values_to="YEARMODA") %>%
  mutate(Year = as.numeric(substr(YEARMODA, 1, 4)),
         Month = as.numeric(substr(YEARMODA, 5, 6)),
         Day = as.numeric(substr(YEARMODA, 7, 8))) %>%
  make_JDay() 
ggplot(data = CKA_Boskoop,
       aes(Pheno_year,
           JDay,
           col = variable)) +
  geom_line() +
  theme_bw(base_size = 15) +
  scale_color_discrete(
    name = "Phenological event",
    labels = c("First bloom",
               "Full bloom",
               "Last bloom")) +
  xlab("Phenological year") +
  ylab("Julian date (day of the year)")

ggplot(data = CKA_Boskoop,
       aes(Pheno_year,
           JDay,
           col = variable)) +
  geom_line() +
  theme_bw(base_size = 15) +
  scale_color_discrete(name = "Phenological event",
                       labels = c("First bloom",
                                  "Full bloom", 
                                  "Last bloom")) +
  xlab("Phenological year") +
  ylab("Julian date (day of the year)") +
  geom_smooth(method = "lm")
## `geom_smooth()` using formula = 'y ~ x'

ggplot(data=CKA_Boskoop,aes(Pheno_year,JDay,col=variable)) +
  geom_smooth() +
  theme_bw(base_size=15) +
  scale_color_discrete(
    name = "Phenological event",
    labels = c("First bloom", "Full bloom", "Last bloom")) +
  xlab("Phenological year") +
  ylab("Julian date (day of the year)") 
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

  1. Evaluate the occurrence of frost events at Klein-Altendorf since 1958. Illustrate this in a plot.
require(Kendall)
Kendall_first <-
  Kendall(x = CKA_Boskoop$Pheno_year[
            which(CKA_Boskoop$variable == "First_bloom")],
          y = CKA_Boskoop$JDay[
            which(CKA_Boskoop$variable == "First_bloom")])

Kendall_full <- 
  Kendall(x = CKA_Boskoop$Pheno_year[
            which(CKA_Boskoop$variable == "Full_bloom")],
          y = CKA_Boskoop$JDay[
            which(CKA_Boskoop$variable == "Full_bloom")])

Kendall_last <- 
  Kendall(x = CKA_Boskoop$Pheno_year[
            which(CKA_Boskoop$variable == "Last_bloom")],
          y = CKA_Boskoop$JDay[
            which(CKA_Boskoop$variable == "Last_bloom")])

Kendall_first
## tau = -0.35, 2-sided pvalue =0.000084584
Kendall_full
## tau = -0.447, 2-sided pvalue =0.00000046775
Kendall_last
## tau = -0.38, 2-sided pvalue =0.000019136

The p-values are all smaller than 0.05, so it is pretty clear that there is a trend.

frost_df = data.frame(
  lower = c(-1000, 0),
  upper = c(0, 1000),
  weight = c(1, 0))

frost_model <- function(x) step_model(x,
                                      frost_df)

hourly <- stack_hourly_temps(CKA_weather,
                             latitude = 50.625)

frost <- tempResponse(hourly,
                      models = c(frost = frost_model))

ggplot(frost,
       aes(End_year,
           frost)) +
  geom_smooth() +
  geom_point() +
  ylim(c(0, NA)) +
  ylab("Frost hours per year") +
  xlab("Year")
## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

  1. Produce an illustration of the relationship between spring frost events and the bloom period of ‘Roter Boskoop’.
frost_model_no_summ <- 
  function(x) step_model(x, 
                         frost_df,
                         summ=FALSE)

hourly$hourtemps[, "frost"] <- frost_model_no_summ(hourly$hourtemps$Temp)

Daily_frost_hours <- aggregate(hourly$hourtemps$frost,
                               by = list(hourly$hourtemps$YEARMODA),
                               FUN = sum)

Daily_frost <- make_JDay(CKA_weather)

Daily_frost[, "Frost_hours"] <- Daily_frost_hours$x

Daily_frost$Frost_hours[which(Daily_frost$Frost_hours == 0)] <- NA

Ribbon_Boskoop <-
  CKA_Boskoop %>%
  select(Pheno_year, variable, JDay) %>%
  pivot_wider(names_from = "variable", values_from = "JDay")

lookup_dates <- Ribbon_Boskoop

row.names(lookup_dates) <- lookup_dates$Pheno_year

Daily_frost[, "First_bloom"]<-
  lookup_dates[as.character(Daily_frost$Year),
               "First_bloom"]

Daily_frost[, "Last_bloom"]<-
  lookup_dates[as.character(Daily_frost$Year),
               "Last_bloom"]

Daily_frost[which(!is.na(Daily_frost$Frost_hours)),
            "Bloom_frost"] <-
  "Before bloom"

Daily_frost[which(Daily_frost$JDay >= Daily_frost$First_bloom),
            "Bloom_frost"]<-
  "During bloom"

Daily_frost[which(Daily_frost$JDay > Daily_frost$Last_bloom),
            "Bloom_frost"]<-
  "After bloom"

Daily_frost[which(Daily_frost$JDay > 180),
            "Bloom_frost"]<-
  "Before bloom"

ggplot(data = Ribbon_Boskoop,
       aes(Pheno_year)) +
  geom_ribbon(aes(ymin = First_bloom, 
                  ymax = Last_bloom),
              fill = "light gray") +
  geom_line(aes(y = Full_bloom)) +
  theme_bw(base_size = 15) +
  xlab("Phenological year") +
  ylab("Julian date (day of the year)") +
  geom_point(data = Daily_frost,
             aes(Year,
                 JDay,
                 size = Frost_hours,
                 col = Bloom_frost),
             alpha = 0.8) + 
  scale_size(range = c(0, 5),
             breaks = c(1, 5, 10, 15, 20),
             labels = c("1", "5", "10", "15", "20"),
             name = "Frost hours") +
  scale_color_manual(
    breaks = c("Before bloom",
               "During bloom",
               "After bloom"),
    values = c("light green",
               "red",
               "light blue"),
    name = "Frost timing") +
  theme_bw(base_size = 15) +
  ylim(c(75, 140))

  1. Evaluate how the risk of spring frost for this cultivar has changed over time. Has there been a significant trend?

    In the viewed period is a trend visible that shows earlier bloom dates, which could increase the risk of frost damage. However, there is no higher greater frost damage risk visible in the data.

Chapter 34 A robust method to estimate future frost risks (?)

There are no tasks for this chapter.

Chapter 35 Major concepts

There are no tasks for this chapter.